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In  this  study,  thermal  spalling  of  high-strength  concrete  structural  elements  has 
been  investigated  by  developing  numerical  models  capable  of  analyzing  the  combined 
effects  of  mass  and  heat  transfer  phenomena  on  thermally  induced  stresses.  Through  the 
use  of  finite  difference  simulations,  this  study  investigates  moisture  effects  on 
thermodynamic  states  of  concrete  at  elevated  temperatures  and  the  influence  of  pore 
pressure  on  development  of  thermal  stress  induced  by  temperature  gradients.  A finite 
element  stress  analysis  model  is  combined  with  finite  difference  simulations  to  predict 
stress  states  capable  of  inducing  spalling  of  the  type  that  has  been  observed  both  in  the 
field  and  in  laboratory  fire  testing  of  concrete  structural  elements.  The  combined  hydro- 
thermo-mechanical  analysis  procedure  presented  here  may  eventually  serve  as  a critical 
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component  of  stress  analysis  used  for  evaluating  the  fire  resistance  of  high-strength 
concrete  structural  systems. 

The  model  of  concrete  exposed  to  fire  that  is  developed  herein  involves  mass  and 
heat  transport  phenomena  in  a multi-phase  continuum.  Simultaneous  flow  of  multiple 
fluid  phases  (i.e.,  liquid  and  gaseous)  is  modeled  using  newly  proposed  constitutive 
relationships  that  are  considered  more  accurate  than  those  previously  available  in  terms 
of  assessment  of  thermodynamic  state  variables  of  concrete  system.  In  addition, 
numerical  simulations  are  used  to  explore  the  effects  that  steel  reinforcing  bars  have  on 
internal  temperature  and  pressure  within  concrete  members. 

Pore  pressure  and  temperature  time  histories  from  hydro-thermal  finite  difference 
simulations  are  presented  and  discussed.  Results  from  the  simulations  yield  an  improved 
understanding  of  the  thermodynamic  state  variables  such  as  pore  pressure,  temperature, 
and  degree  of  liquid  water  saturation.  Selection  of  constitutive  models  to  describe  the 
flow  characteristics  of  concrete  and  the  procedures  implemented  for  the  creation  of  the 
concrete  models  are  also  discussed.  Modeling  aspects  such  as  meshing  techniques, 
selection  of  initial  conditions,  and  definition  of  boundary  conditions  are  discussed.  The 
modeling  techniques  used  to  represent  radiant  heat  flux  boundary  conditions  and  volume- 
averaging processes  are  explained  in  detail. 

Finally,  a theory  of  stress  superposition  is  developed  and  subsequently  used  to 
evaluate  stress  states  so  as  to  determine  the  controlling  contributor  to  thermal  spalling  of 
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high-strength  concrete,  e.g.,  hydrostatic  pore-pressures  and  thermal  gradient  stresses, 
under  extreme  thermal  loading  conditions. 
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CHAPTER  1 
INTRODUCTION 


The  thermo-mechanical  behavior  of  concrete  structures  exposed  to  high  temperature 
environments  has  become  very  important  for  many  engineering  problems  in  recent  years. 
Assessment  of  the  fire  resistance  of  buildings,  bridges  and  tunnels,  and  assessment  of  fire 
damage  in  structures  represent  a few  example  situations.  In  these  areas,  many  engineers 
consider  understanding  the  fire  resistant  characteristics  of  concrete  to  be  a critical 
component  of  accurate  structural  analysis  and  thus  safe  yet  economical  designs. 

Economically,  concrete  has  better  fire-resistant  properties  than  other  building 
materials  and  thus,  it  provides  effective  shielding  to  materials  such  as  steel  against 
serious  damage  during  fires.  However,  in  an  extremely  rapid  heating  condition  resulting 
from  an  elevated  temperature  environment,  sudden  explosive  material  failure  of  concrete, 
i.e.,  spalling,  may  occur,  rendering  it  incapable  of  protecting  steel  reinforcement  and 
structurally  unable  to  carry  load. 

If  a reinforced  concrete  column  in  an  extreme  fire  environment  undergoes  spalling 
at  its  surface,  then  the  steel  reinforcement  can  be  directly  exposed  to  fire.  Strength  of  the 
steel  at  higher  temperatures  will  be  progressively  reduced.  Sudden  loss  of  gross  sectional 
area  in  the  concrete  column  due  to  spalling  can  also  result  in  a reduction  of  column 
stiffness.  As  spalling  progresses  into  deeper  regions  of  the  structural  member’s  core,  the 
concrete  column  may  become  incapable  of  carrying  the  acting  loads  as  severe  loss  of 
strength  occurs. 
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The  protection  of  steel  beams  and  girders  by  a concrete  slab  is  also  important  in 
assessing  the  fire  safety  of  structures.  With  a fire  above  the  slab,  the  structural 
strength  of  beams  can  be  significantly  changed  by  the  fire  proofing  capability  of  the 
slab  that  provides  shielding  of  the  beams  against  heat. 

Fire  in  a tunnel  may  cause  serious  safety  problems  of  structural  integrity  and 
thus,  economic  loss  and  loss  of  life.  One  example  case  involved  a fire  in  the 
Eurotunnel.  During  a fire  incidence  in  the  channel  tunnel  (Ulm  et  al.  1999a),  severe 
spalling  occurred  in  the  tunnel  lining.  Almost  200  mm  deep  surface  was  delaminated 
and  the  steel  reinforcement  lost  its  40%  yield  strength  due  to  heat  reaching  1 100  °C. 
Concerned  with  the  possibility  of  large  scale  structural  damage,  the  necessity  of 
understanding  a potential  structural  failure  mechanism  in  combination  with  the 
material  behavior  of  concrete  at  high  temperatures  has  been  emphasized  in  the  study 
of  tunnel  fire  incidents. 

Thus,  the  focus  of  this  dissertation  is  the  development  of  more  reliable  solution 
methods  to  study  thermo-mechanical  behavior  of  concrete.  The  development  of  an 
analytical  model  including  the  formulation  and  the  application  of  the  finite  element 
method  to  calculate  stress  distribution  is  developed  so  that  the  composition  of  stress 
produced  in  concrete  structures  exposed  to  fire  may  be  determined.  Once  such  a 
model  has  been  developed  and  validated  for  a reasonable  range  of  applications,  it 
would  be  used  to  determine  fire  safety  assessment  in  existing  buildings  and  develop 
design  guidelines  and  specifications  for  the  use  of  concrete  in  regard  to  fire  safety. 
Therefore,  a numerical  tool  including  material  models  capable  of  representing  the 
relevant  physical  phenomena  that  occur  in  heated  concrete  needs  to  be  developed. 
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1.1  Material  Response  of  Concrete  to  Heating 

Concrete  is  physically  a three-phase  material  as  shown  in  Figure  1-1.  The  first 
phase  considered  in  this  study  is  the  solid  phase.  It  consists  of  the  cementitious  matrix 
made  of  aggregates,  hydrated  cement,  and  a bonding  interface  between  the  hydrated 
cement  and  the  aggregates.  The  second  phase  includes  the  liquid  phase  either  in  the 
form  of  adsorbed  water  within  the  solid  phase  (e.g.  chemically  bound  water  in  the 
hydrated  cement)  or  in  the  form  of  pore  water  (e.g.  free  water  or  evaporable  water  in 
the  pores).  The  gaseous  phase  (i.e.,  the  mixture  of  dry  air  and  water  vapor)  occupying 
the  pores  makes  up  the  third  phase. 


Gaseous  phase  (air- 

Liquid  phase  (e.g.  liquid 
Entrapped  gas 

Figure  1-1.  A schematic  drawing  of  constituents  of  concrete 

When  concrete  is  exposed  to  fire,  capillary  water  content  evaporates.  The 
generated  water  vapor  migrates  along  pressure  gradients  within  the  concrete  matrix 
toward  both  the  surface  of  the  matrix  and  deeper  into  the  matrix.  As  the  temperature 
keeps  rising  on  the  exposed  surface,  desaturation  inside  the  concrete  intensifies  with 
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high  heat  flow,  and  then,  a major  portion  of  the  released  vapor  is  driven  by  a large 
pressure  gradient  toward  cooler  regions.  The  vapor  condenses  in  the  pores  of  the 
cooler  regions  as  the  pressure  buildup  levels  off.  The  coupled  moisture  and  heat  flow 
resistance  along  the  flow  path  increases  (i.e.,  permeability  decreases)  and  the 
temperature  redistribution  inside  the  concrete  is  affected  by  the  flow  (i.e., 
macroscopic  conductivity  changes).  The  multiphase  transport  and  phase  change 
mechanisms  of  the  water  content  generate  a coupled  heat  and  mass  transfer 
phenomenon  in  the  porous  medium. 

Despite  the  significant  involvement  of  mass  transport  phenomena  in  the 
temperature  distributions  that  arise,  the  multi-phase  characteristics  of  concrete  have 
not  been  emphasized  in  the  literature  dealing  with  heat  transfer  analysis  in  concrete. 
As  a result,  the  coupled  heat  and  mass  transport  has  generally  been  neglected  or 
simply  assumed  to  be  negligible  in  assessing  the  temperature  gradients  and  thus 
computation  of  thermal  gradient  stresses. 

In  this  study,  the  moisture  contents  and  their  transport  mechanisms  occurring  in 
heated  concrete  are  studied  and  the  coupling  effects  are  considered.  Thus,  thermal 
gradients  stresses,  which  are  suspected  as  one  of  the  major  causes  of  eventual 
material  failure,  are  quantified  based  on  the  thermal  gradients  computed  by  coupled 
moisture  and  heat  transfer. 

1.2  Failure  Mechanisms  of  Concrete  at  Elevated  Temperatures 

When  concrete  structures  are  exposed  to  high  temperature  sources  such  as  fire, 
material  failure  associated  with  rapid  heating  can  occur.  The  term  “spalling”  is  used 
to  describe  the  separation  of  concrete  from  the  surface  of  a concrete  element — such  as 
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columns,  beams,  slabs  and  walls — when  concrete  structural  elements  are  rapidly 
heated.  This  material  failure  can  seriously  affect  fire  resistance  of  a structure  and  thus 
structural  performance.  The  separation  process  can  be  classified  into  the  following 
categories: 

Spalling.  This  is  defined  as  a violent  form  of  material  failure  that  occurs  at  an 
early  stage  of  heating  and  may  result  in  extensive  damage  or  complete  failure  of  the 
concrete  element. 

Sloughing  off.  This  is  defined  as  a progressive  form  of  breakdown,  which 
involves  gradual  separation  over  longer  periods  of  time  of  material  from  the  concrete 
element.  It  may  continue  slowly  through  the  later  stage  of  heating. 

The  severity  and  sudden  failure  mechanism  of  spalling  is  regarded  as  being  more 
serious  in  the  assessment  of  structural  safety  and  thus  requires  additional  attention. 

1.3  Two  Hypotheses  in  Analysis  of  Spalling 

The  causes  of  spalling  have  been  studied  for  several  decades.  Two  hypotheses 
have  been  proposed  to  explain  explosive  spalling.  The  first  hypothesis  is  that  heating 
produces  very  high  water  vapor  pressure  in  the  concrete  pores  (Harmathy  1964), 
which  results  in  large  tensile  stresses  and  consequent  tensile  failure.  In  many 
numerical  approaches  based  on  this  hypothesis  (Sahota  and  Pagni  1979,  Kodres  1996, 
Tenchev  et  al.  2001),  a microscopic  pore  pressure  buildup  was  treated  as  a 
macroscopic  stress  component  of  the  heated  concrete.  With  a macroscopically 
averaging  process,  contribution  of  the  microscopic  thermodynamic  pressure  can  be 
substantially  reduced  to  an  actual  stress  state  that  should  be  used  in  the  determination 
of  potential  macroscopic  material  failure.  It  has  previously  been  unknown  as  to 
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whether  pore  pressure  is  a major  contributor  to  spalling.  The  key  material  parameters 
such  as  permeability  and  conductivity  of  concrete,  which  determine  flows  of  the  fluid 
phases  and  which  can  be  determining  factors  in  the  coupling  effects  between  heat  and 
moisture  flows,  have  to  be  carefully  chosen  in  computational  analyses. 

The  second  is  that  sudden  localized  heating  causes  thermal  expansion,  which  is 
restrained  by  the  adjacent  cool  concrete.  This  phenomenon  produces  a large 
compressive  stress,  which  is  suspected  to  cause  crushing  of  the  concrete  and  buckling 
of  the  surface  layer.  Although  high  pore  pressure  develops  inside  concrete,  it  only 
serves  to  create  cracks.  Nevertheless,  once  crack  opening  starts,  the  pore  pressure  is 
assumed  to  be  dissipated  and  drop  nearly  to  zero.  Furthermore,  on  the  basis  of 
fracture  mechanics,  brittle  fracture  due  to  sudden  release  of  the  potential  energy  of 
the  thermal  stress  has  been  explained  as  a major  cause  to  explosive  spalling  (Bazant 
1997)  where  concrete  flaws  may  act  as  fracture  nuclei.  The  medium  was  assumed  to 
be  a homogeneous  continuum  in  the  sense  that  the  size  of  a dominant  flaw  is  large  in 
comparison  with  the  characteristic  microstructural  dimensions  of  the  material 
constituents.  Initiation  and  the  process  of  material  failure  of  concrete  then  occur  in  the 
influence  of  applied  loads,  flaw  geometry,  size  of  the  structural  element  and  material 
behavior  in  crack  propagation  from  the  dominant  flaw  in  concrete. 

One  thing,  however,  has  been  unclear  as  to  whether  thermal  expansion  is  more 
dominant  a mechanism  than  thermal  gradient  effects  when  coupling  effects  between 
heat  flow  and  steam  induced  by  evaporation  of  the  free  water  affect  the  temperature 
distribution.  Temperature  rise  can  be  significantly  different  if  fluid  phases  (i.e.,  liquid 
and  gaseous  phases)  in  concrete  pores  are  included  in  the  analysis  of  heat  transfer  The 
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role  of  the  fluid  phases  in  development  of  the  stress  field,  which  is  possibly 
interrelated  with  crack  propagation  at  elevated  temperatures,  needs  to  be  included. 
Thus,  influence  of  pore  pressure  and  its  effects  in  cracking  of  heated  concrete  should 
be  considered  for  the  relevant  thermodynamic  state  in  heated  concrete.  In  this 
dissertation,  pore  pressure  is  examined  to  determine  whether  it  is  a significant 
contributor  to  the  stress  state  and  to  material  failure. 

1.4  Moisture  Clog  Spalling 

Spalling  has  generally  been  found  to  be  more  severe  (as  observed  in  actual  fires 
and  during  laboratory  testing  of  specimens)  in  high  strength  concrete  (HSC)  than  in 
normal  strength  concrete  (NSC)  (Sullivan  2001).  To  investigate  a process  known  as 
moisture  clog  spalling  (Harmathy  1964),  flow  of  both  fluid  phases  (i.e.,  liquid  and 
gaseous  phases)  needs  to  be  studied.  The  process  of  vaporization  and  condensation 
continues  during  heating  and  eventually  a “moisture  clog”,  i.e.  a region  of 
significantly  increased  pore  liquid  saturation,  has  the  effect  of  significantly  reducing 
the  permeability  of  the  porous  matrix  to  gas  (steam)  flow.  Consequently,  the  process 
results  in  a sharp  rise  in  pore  pressure  at  that  location,  and  induces  large  differentials 
of  the  temperature  field  of  the  matrix,  which  consequentially  causes  a significant 
change  of  the  development  of  temperature  gradients.  If  the  resistance  of  the  pore 
network  to  moisture  flow  is  overcome  by  the  pressure  gradient,  the  condensed 
moisture  layer  begins  to  move  deeper. 

In  this  study,  attention  is  given  to  this  phenomenon  that  generally  has  been 
suspected  to  be  responsible  for  severe  spalling  in  HSC.  Because  of  the  lower 
permeability  of  HSC,  pore  pressure  induced  by  water  vaporization  is  unable  to 


8 


dissipate  quickly.  Several  factors  that  affect  moisture  flow  including  the 
permeability,  conductivity,  porosity,  and  initial  liquid  saturation  level  of  concrete  will 
be  studied  in  this  dissertation.  A more  detailed  explanation  and  prediction  of  pore 
pressure  buildup  will  be  addressed  by  investigating  multi-phase  mass  flow  that 
nonlinearly  couples  with  heat  transfer  mechanisms  such  as  internal  convection  and 
conduction. 


1.5  Overview  of  the  Dissertation 

The  need  for  more  reliable  solution  methods  to  study  the  behavior  of  the 
concrete  structures  exposed  to  fire  has  prompted  the  development  of  a coupled  multi- 
phase analytical  model  in  the  present  study.  Capable  of  analyzing  the  thermally 
induced  stress  state  under  consideration  of  heat  and  mass  transfer,  a computational 
model  forms  the  basis  for  this  research. 

As  a parametric  study,  the  first  part  of  this  study  focuses  on  the  use  and 
modification  of  an  existing  finite  difference  based  heat  and  mass  transport  simulation 
code  in  order  to  examine  pore  pressures  quantitatively  and  determine  temperature 
distributions  in  concrete  exposed  to  fire.  In  the  modeling  of  moisture  movement 
through  concrete,  most  important  among  the  issues  are  the  development  of  a new 
relative  permeability  function  for  concrete  (presented  herein)  and  consideration  of  the 
slippage  effect  in  gas  flow  through  heated  concrete.  Because  heat  and  mass  transfer  in 
concrete  has  significant  dependence  on  the  relevant  material  characteristics, 
particularly  on  permeability,  the  phenomena  such  as  relative  permeability  and  gas  slip 
flow  effects  need  to  be  considered  in  determination  of  thermodynamic  state  variables 
such  as  pore  pressure  and  temperature.  In  consideration  of  the  fundamental 
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relationships  among  these  parameters,  a numerical  solution  is  sought  and  compared  to 
experimental  tests.  Using  a finite  difference  analysis  (FDA)  model,  key  field 
parameters  for  the  modeling  are  identified. 

In  the  second  part  of  the  dissertation,  multi-dimensional  finite  element  stress 
analysis  models  are  developed,  which  utilize  the  results  obtained  from  the  finite 
difference  analysis.  Using  the  theory  of  thermodynamics,  practical  assumptions  used 
in  the  study  are  discussed  and  are  validated.  Prediction  of  a combined  stress  state 
developed  in  the  heated  concrete  structural  members  requires  simulation  of  the 
structural  boundary  and  thermal  loading  conditions.  Thus,  modeling  techniques  are 
developed  that  take  into  account  realistic  structural  boundary  and  initial  conditions. 

A mathematical  model,  which  can  be  applied  to  a wide  variety  of  three-dimensional 
configurations  of  concrete  structures  such  as  reinforced  concrete  walls,  slabs, 
columns  and  beams  under  different  boundary  conditions,  will  be  presented  in  detail. 
Further  mathematical  discussion  is  made  regarding  energy  balance  in  a non- 
isothermal  partially  saturated  porous  medium  as  a starting  point  for  future  research 
activities. 

Finally,  thermal  spalling  of  structural  concrete  elements  is  studied  and  analyzed 
by  making  use  of  the  newly  proposed  models.  Effects  of  the  coupled  heat  and  mass 
transfer  mechanism  on  the  prediction  of  combined  stress  states  are  studied 
quantitatively.  With  a better  knowledge  of  the  thermodynamic  state,  thermally 
induced  stresses  are  then  investigated  using  the  theory  of  elasticity.  Therefore,  the 
contribution  of  the  pore  pressure  buildup  to  spalling  can  be  identified.  In  addition,  the 
rate  of  penetration  of  the  “moisture  clog”  and  the  locations  of  thermodynamic  pore- 
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pressure  peaks  under  different  thermal  loading  conditions  are  discussed,  which  can  be 
crucial  to  understanding  material  failure  mechanisms  such  as  spalling. 


CHAPTER  2 
LITERATURE  REVIEW 


For  many  years,  heat  and  mass  transfer  in  heated  concrete  has  been  investigated 
both  theoretically  and  experimentally  as  numerous  engineers  and  scientists  have  studied 
material  behaviors  of  porous  media  at  elevated  temperatures.  Generally,  three  different 
approaches  to  the  problem  can  be  observed  (Lewis  and  Schrefler  2000): 

Microscopic  level.  “The  scale  of  nonhomogeneity  is  of  the  same  order  of 
magnitude  as  the  dimensions  of  a pore  or  a grain.  Therefore,  a mathematical  point  within 
a single  phase  and  the  field  variables  describing  the  status  of  a phase  are  defined  only  at 
the  points  occupied  by  that  phase”  (pp.  10). 

Macroscopic  level.  “The  real  multiphase  system  that  occupies  the  porous  medium 
domain  is  replaced  by  a model  at  the  level  of  interest  of  continuum  mechanics  where  a 
continuous  distribution  of  the  constituents  through  a macroscopic  control  space  is  treated 
using  an  averaging  process”  (pp.  1 0). 

Megascopic  level.  “Selected  macroscopic  nonhomogeneities  are  eliminated  by 
averaging  and/or  by  using  a mathematical  model  that  is  stated  in  a domain,  which  has 
fewer  dimensions  than  the  real  domain  (e.g.,  a two  dimensional  problem  with  field-values 
averaged  over  the  thickness)”  (pp.10). 

In  theory,  the  governing  equations  that  describe  heat  and  mass  transport  phenomena 
in  concrete  can  be  written  at  the  microscopic  level  where  only  a mathematical  point 
within  a considered  phase  present  in  the  domain  is  focused  on.  However,  the  surface 
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boundaries  of  the  phases  in  this  type  of  system  description  are  too  complex  to  be 
described.  The  pore  structure  of  concrete  has  an  extremely  complicated  geometry  so 
that  a mechanical  description  of  the  problem  at  the  microscopic  level  is  difficult.  The 
state  variables  (e.g.,  velocity,  pore-pressure  in  the  fluid  phase,  temperature  and  stress 
in  the  solid  skeleton)  within  a phase  vary  from  point  to  point. 

For  engineering  purposes,  a model  at  the  macroscopic  level  or  above  is  of 
interest,  where  measurable,  continuous  and  differentiable  extensive  quantities 
(quantities  that  can  be  additive  over  volumes  with  mass,  momentum  and  energy)  can 
be  determined  and  boundary  value  problems  can  be  stated  and  solved.  Such  models 
describe  interacting  constituents  such  as  the  solid,  liquid  and  gaseous  phases  that 
occupy  the  entire  control  space  of  the  porous  medium.  The  description  of  concrete  as 
a multiphase  system  made  of  those  continuous  bodies  is  based  on  averaging  theories 
of  field  variables  that  dictate  the  state  of  the  system.  By  using  the  averaging  theory 
based  on  spatial  averaging  processes,  macroscopic  variables  that  correspond  to  real 
laboratory  measurable  quantities  are  used  for  developing  the  governing  equations. 

The  main  goals  of  this  chapter  are  to  review  the  continuum  approaches  used  in 
the  past  and  to  lead  the  readers  to  the  macroscopic  level  of  describing  phenomena  in 
mass  and  heat  transfer  in  porous  media.  The  governing  equations  of  a partially 
saturated  concrete  porous  medium  at  elevated  temperature  are  then  derived  and 
discussed  in  the  following  chapter. 

2.1  Relevant  Flow  Properties 

In  modeling  simultaneous  multiphase  fluid  flows  in  heated  concrete,  differences 
between  permeability  to  gaseous  phase  and  permeability  to  liquid  phase  of  the 
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concrete  have  been  overlooked  in  many  past  numerical  studies  (Huang  et  al.  1979, 
Sahota  and  Pagni  1979,  Hurst  and  Ahmed  1998).  Significant  changes  in  gas 
permeability  of  concrete — when  moisture  is  vaporized  and  redistributed  in  elevated 
temperatures — were  not  considered  in  the  numerical  predictions  of  pore  pressure 
buildup.  Furthermore,  gas  permeability  of  concrete  was  arbitrarily  chosen  in  a few 
cases  of  numerical  simulations. 

Relative  permeability  is  a concept  that  is  very  important  for  the  simulations 
considered  in  this  study.  It  includes  macroscopic  mathematical  measurement 
describing  the  multiphase  flow  of  fluids  in  porous  media  (Corey  et  al.  1956,  Burdine 
1953)  and  represents  the  fluid-fluid  interactions  during  multiphase  flow  (Li  and 
Home  2001).  The  relative  permeability  of  a porous  network  can  be  thought  of  as  a 
fraction  of  the  intrinsic  permeability  for  the  phase  of  interest.  When  two  or  more 
phases  flow  simultaneously  inside  concrete,  each  phase  will  develop  its  own  flow 
network  to  which  a permeability  can  be  assigned  in  much  the  same  way  as  for  single- 
phase flow. 

A new  relative  permeability  function  for  high  strength  concrete  (concrete  that 
has  a low  water-cement  ratio  and  contains  silica  fume  (or  fly  ash)  and 
superplasticizer)  as  well  as  normal  strength  concrete  (concrete  that  has  a relatively 
higher  water-cement  ratio  and  a compressive  strength,  fc , of  6000  psi  or  below)  is 
developed  and  then  implemented  as  a new  function  into  the  TOUGH2  code  (Pruess 
1991).  Reflecting  the  manner  in  which  flowing  fluid  phases  are  distributed  within 
concrete,  the  idea  of  relative  permeability  was  adapted  to  account  for  the  wetting 
state,  pore  structure  and  saturation  history  of  concrete  structures  during  mass  and  heat 
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transfer.  Due  to  moisture  clog  formation  during  mass  and  heat  transfer,  flow  of  the 
gaseous  phase  is  affected  and  inhibited.  Relative  permeability  functions  written  as 
functions  of  degree  of  phase  saturation  were  thus  developed.  This  is  a new 
application  for  hysteretic  constitutive  relations  governing  multiphase  flow  inside 
concrete  exposed  to  rapid  heating.  In  a recent  literature  search,  the  concept  of  relative 
permeability  was  verified  by  experimental  measurements  of  gas  permeability  of 
partially  saturated  concrete  (Jacobs  1998).  Chung  and  Consolazio  (2003)  showed  that 
numerical  computation  of  mass  flow  rate  in  consideration  of  fluids’  relative 
interaction  inside  heated  concrete  was  significantly  changed  by  the  relative 
permeability  implementation.  Changes  in  prediction  of  mass  transport  and 
temperature  gradient  development  indicated  the  importance  of  the  relative 
permeability  concept  in  modeling  coupled  mass  and  heat  transfer  phenomena  in 
heated  concrete.  In  addition  to  phase  interference  modeling,  gas  phase  slip-flow 
effects  are  also  important  in  the  present  context. 

The  slip  flow  phenomenon  that  arises  during  gas  movement  in  concrete  can  be 
quite  substantial  for  low  permeability  materials  (Klinkenberg  1941)  such  as  high 
strength  concrete  (HSC).  Gas  permeability  of  concrete  may  vary  approximately  two 
orders  of  magnitude  due  to  the  slip  effects.  Thus,  when  interpreting  data  from  gas 
permeability  testing  of  concrete,  one  must  be  able  to  isolate  and  “remove”  the 
Klinkenberg  (slip  flow)  effects  in  order  to  determine  the  fundamental  intrinsic 
permeability  of  the  concrete. 

Experimental  studies  by  McVay  and  Rish  (1995)  revealed  that  slip  flow,  not 
laminar  flow,  occurs  when  superheated  steam  and  nitrogen  gas  are  passed  through 
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cement  mortar.  They  demonstrated  a procedure  involving  testing  at  various  levels  of 
gas  pressure  that  quantifies  the  slip  flow  effect  with  b (Klinkenberg’s  constant)  data 
and  evaluates  Kg  (the  fundamental  intrinsic  permeability)  data  for  the  tested 

material.  For  numerical  modeling  purposes,  the  parameters  Kg  and  b are  used  for 

gas  flow  calculations  and  the  liquid  flow  parameter  K(  (the  water  permeability)  is 
used  to  compute  liquid  flow  rates  (Bamforth  1987a). 

However,  Kf  has  often  been  treated  as  being  equivalent  to  Kg  , the  fundamental 

intrinsic  permeability  of  the  concrete  (Bamforth  1987b).  This  results  in  an 
overestimation  of  the  slip  flow  phenomenon  in  gas  permeability  of  concrete.  Dhir  et 
al.  (1989)  conducted  experimental  tests  to  identify  the  relationship  between  air  and 
water  permeability  of  concrete.  They  demonstrated  that  even  when  the  gas  slippage 
and  temperature  effects  are  considered  in  the  calculation  of  the  fundamental  intrinsic 
permeability  of  concrete,  the  derived  intrinsic  permeability  of  concrete  could  still 
vary  with  respect  to  whether  air  or  water  was  used  in  the  measurement  process.  The 
difference  in  experimental  results  from  gas  and  liquid  permeability  of  concrete  was 
also  reported  by  Whiting  (1988).  The  departure  from  the  concept  of  a single  absolute 
value  of  concrete  permeability  is  believed  to  be  caused  by  physical  (e.g.,  dilation  of 
cement  hydrates  in  contact  with  water)  and  chemical  (e.g.,  hydration  of  cement  paste) 
interactions  during  the  long  testing  periods  used  in  water  permeability  testing  of 
concrete. 

Chung  and  Consolazio  (2003)  developed  an  interpretation  method  for  using 
experimental  data  from  gas  permeability  testing  of  concrete  to  determine  the 
Klinkenberg  constant  b that  quantifies  just  how  sensitive  the  “apparent”  gas 
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permeability  of  the  material  is  to  pressure  changes.  Emphasis  was  on  properly 
accounting  for  slip  flow  in  numerical  modeling  of  mass  transport  inside  heated 
concrete,  particularly  HSC  (because  its  lower  permeability  is  more  affected  by 
pressure  changes  than  normal  strength  concrete  (NSC)). 

2.2  Material  Characteristics  of  High  Strength  vs.  Normal  Strength  Concrete 

Pore  structural  network  characteristics  of  concrete  have  generally  been 
suspected  to  result  in  the  risk  of  explosive  thermal  spalling,  especially  for  structural 
members  made  of  high  performance  or  strength  concrete.  For  example,  permeability 
in  HSC  is  more  severely  affected  by  changes  of  the  degree  of  saturation  as  the 
evaporation-condensation  process  forms  moisture  clog. 

Huang  et  al.  (1979)  pointed  out  that  the  distribution  of  pore  size  might  be 
considered  the  most  important  parameter  affecting  moisture  transfer  in  a porous 
medium.  The  characteristics  of  concrete  related  to  the  distribution  of  pore  size  were 
identified  later  by  Baroghel-Bouny  and  Chaussadent  (1996).  They  experimentally 
verified  that  classical  pore  structure  investigation  methods  such  as  mercury  intrusion 
or  nitrogen  adsorption  were  not  able  to  provide  accurate  pore  network  characteristics 
of  cement-based  materials,  particularly  for  high-performance  concrete. 

Additional  efforts  to  understand  pore  structures  have  focused  on  using  sorption 
curves  to  calculate  the  local  equilibrium  moisture  content  (Hansen  1989).  Water 
vapor  desorption  and  adsorption  isotherm  tests  have  been  conducted  for  both  NSC 
and  HSC.  The  experimental  water  vapor  sorption  isotherms  (hydrometrical 
equilibrium  at  a constant  temperature)  showed  that  more  severe  material  hysteresis  in 
isotherms  was  observed  in  HSC.  Kjellsen  and  Atlassi  (1999)  also  investigated 
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influence  of  silica  fume  on  the  pore  structure  of  mortar,  using  water  vapor  sorption 
isotherms.  It  was  concluded  that  the  characteristics  of  pore  size  distribution  and  the 
ratio  of  pore  volume  to  total  volume  of  the  concrete  specimen  were  the  key  to  cause 
such  a significant  hysteresis  in  isotherms  of  HSC.  The  test  results  proved  that  the 
initial  degree  of  saturation  at  ambient  temperatures  and  the  structural  pore-network 
(i.e.,  permeability  and  porosity)  of  HSC  are  different  from  NSC. 

These  differences  are  thought  to  be  related  to  the  fact  that  behavior  of  saturated 
HSC  subjected  to  intensive  heating  has  been  reported  as  prone  to  explosive  spalling 
(Sanjayan  and  Stocks  1993,  Chan  et  al.  1999).  Chan  et  al.  (1999)  also  showed  that 
moisture  content  had  a dominant  influence  on  spalling,  which  was  concluded  to  be 
related  to  vapor  pressure  buildup.  The  dependence  of  spalling  on  moisture  content 
and  strength  of  concrete  at  elevated  temperatures  was  confirmed.  However,  the 
effects  of  the  two  factors  on  thermally  induced  stress,  and  thus,  material  failure  in  a 
concrete  system  remained  unexplained. 

2.3  Material  Properties  of  Heated  Concrete 

Generally,  strength  of  concrete  at  elevated  temperatures  is  known  to  diminish. 
Many  studies  have  been  conducted  to  determine  thermal  and  mechanical  properties  of 
concrete  at  elevated  temperatures  for  use  in  fire  resistance  calculations  (Castillo  and 
Durrani  1999,  Lie  and  Kodur  1996,  Phan  and  Carino  1998,  Reis  et  al.  2001,  and 
Janotka  and  Bagel  2002). 

The  influence  of  temperature  on  the  mechanical  properties  such  as  compressive 
strength  and  modulus  of  elasticity  was  claimed  to  be  greater  than  the  influence  on 
thermal  properties  of  thermal  conductivity,  specific  heat  and  thermal  expansion  (Lie 
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and  Kodur  1996).  Phan  and  Carino  (1998)  summarized  experimentally  tested  elastic 
moduli  and  compressive  strengths  of  HSC  and  NSC  exposed  to  rapid  heating.  In  their 
report,  HSC  and  NSC  exhibited  differences  in  the  mechanical  properties  in  the 
temperature  range  between  room  temperature  and  450  °C  and  a higher  potential  of 
HSC  to  fail  by  explosive  spalling  was  noted.  Reis  et  al.  (2001)  endeavored  to  measure 
the  reduction  of  fire  resistance  of  HSC.  HSC  was  found  to  be  more  severely  affected 
in  its  mechanical  behavior  when  exposed  to  fire.  For  a remedy  to  explosive  spalling, 
they  suggested  that  using  steel  fiber  reinforcement  could  improve  the  performance  of 
HSC  and  NSC  elements  in  fire.  Janotka  and  Bagel  (2002)  conducted  experimental 
tests  of  thermo-mechanical  properties  of  HSC  at  temperatures  up  to  800  °C.  The  pore 
structure  of  HSC  was  found  to  change  significantly  and  resulted  in  strength  reduction 
and  permeability  increase. 

2.4  Heat  and  Mass  Transport  Phenomena  in  Heated  Concrete 

In  a macroscopic  viewpoint,  the  theory  of  irreversible  process  ( the  second  law 
of  thermodynamics)  can  be  used  to  identify  the  driving  forces,  based  upon  relating  the 
entropy  of  the  porous  system  to  the  heat  transfer  rate,  the  rate  of  mass  flow  and  the 
distribution  of  temperature  and  pore  pressure  (Luikov  1966). 

Elevated  temperature  in  the  porous  system  will  induce  a temperature  gradient 
that  causes  a pressure  gradient  (Sahota  and  Pagni  1979,  Huang  1979,  Wang  and  Yu. 
1988).  Sahota  and  Pagni  (1979)  obtained  a quasi-static  solution  of  a two-dimensional 
transient  conduction  problem  with  mixed  boundary  conditions,  assuming  that 
property  variations  and  moisture  transfer  were  negligible.  Huang  (1979)  investigated 
the  effects  of  mass  concentration,  temperature  and  pore  pressure  gradients  on  energy 
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and  mass  transfer.  Moisture  migration  through  porous  media  was  assumed  to 
predominantly  occur  in  the  gaseous  phase  driven  by  pressure  and  mass  concentration 
gradients.  Liquid  phase  movement  was  neglected  because  of  its  very  small  Darcy’s 
coefficient  ( coefficient  of  permeability,  also  called  hydraulic  conductivity ) compared 
to  that  of  gases.  The  evaporation-condensation  mechanism  was  concluded  to  be  a 
governing  factor  in  the  drying  process  (Huang  et  al.  1979). 

Abdel-Rahman  and  Ahmed  (1996)  developed  a numerical  solution  method  in 
order  to  compute  the  pore  pressure  fluctuations  and  the  influence  of  the  pore  pressure 
on  spalling  of  concrete  at  elevated  temperatures.  Kodres  (1996)  hypothesized  that 
high  pore  pressures  and  pressure  gradients  developed  by  vaporized  water  content  and 
from  air  being  compressed  by  liquid  water  could  be  sufficient  to  cause  a tensile 
failure  in  concrete  when  water  and  air  were  migrating  through  a heated  section  of 
concrete.  Thermal  expansion  of  liquid  pore  water  and  transport  of  water  were  also 
found  to  affect  the  build-up  of  pore  gas  pressure  by  Kalifa  et  al.  (2000).  Their 
experiments  demonstrated  a correlation  between  pressure  peaks  and  temperature 
distributions  (e.g.,  a plateau  in  temperature  curves  was  found  where  the  peak  pore 
pressures  developed).  The  two  hypotheses  confirmed  that  (1)  a drying  front  is 
preceded  by  a saturated  layer  acting  as  moisture  clog  and  (2)  flow  of  fluids  in  pores 
affects  temperature  re-distribution  inside  the  system 

2.5  Stress  analysis  of  concrete  exposed  to  fire 

For  the  fire  incident  that  occurred  on  November  1 8,  1996  in  the  Channel  Tunnel 
connecting  England  and  France,  Ulm  et  al.  (1999a  and  b)  studied  thermal  spalling  to 
evaluate  the  total  stress  state  that  might  have  developed  during  the  fire.  In  their 


20 


hypothesis,  the  magnitude  of  the  pore  pressure  buildup  may  not  be  large  enough  to 
cause  tensile  failure  of  concrete.  It  was  believed  that  the  pore  pressure  would  be 
quickly  dissipated  (or  reduced)  with  crack  opening  that  would  prevent  pressure  build 
up.  As  a result,  the  pore  pressure  buildup  was  claimed  as  a secondary  factor  in 
determination  of  the  spalling  mechanism  with  expansion  of  a growing  crack. 

Moisture  movement  was  not  considered  in  the  heat  transfer  analysis.  A finite  element 
model  (Ulm  et  al.  1999a)  was  developed  with  their  newly  proposed  constitutive 
model,  such  as  isotropic  chemoplastic  softening  model  (Ulm  at  al.  1999b),  for  a 
heated  concrete.  Their  study  in  brittle  fracture  due  to  restrained  thermal  dilatation 
concluded  that  material  changes  caused  by  chemical  interactions  (e.g.,  dehydration) 
in  the  constituents  of  concrete  were  more  influential  than  a coupling  effect  between 
heat  transfer  and  fluid  transport  in  the  thermal  stress  analysis.  Temperature  re- 
distributions induced  by  the  effects  of  multi-phase  interactions  of  heat  transfer  were 
not  considered. 

Baggio  et  al.  (1995)  conducted  a computational  analysis  of  mechanical  behavior 
of  concrete  structures  combined  with  thermo-hygrometric  effects  of  the  drying 
process  in  ambient  temperatures.  They  emphasized  that  a thorough  understanding  of 
the  heat  and  moisture  transfer  phenomena  occurring  inside  heated  concrete  was 
required  to  more  accurately  evaluate  thermal  and  moisture  transmission  behavior  and 
the  related  stresses  and  strains  of  the  porous  building  material.  With  an  understanding 
of  the  influence  of  the  fluid  transport  mechanisms  on  heat  transfer  phenomena,  the 
efforts  are  extended  in  the  present  study  to  modeling  the  coupled  effects  and 
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analyzing  a combined  stress  state  of  the  heated  concrete,  which  results  from  both 
thermo-mechanical  stresses  and  pore  pressure. 


CHAPTER  3 

MACROSCOPIC  DESCRIPTION  OF  TRANSPORT  PHENOMENA  IN  HEATED 

CONCRETE 

The  objective  of  this  chapter  is  to  develop  a mathematical  model  of  a multi-phase 
continuum  that  describes  heat  transfer  and  mass  transport  phenomena  in  heated  concrete 
at  the  macroscopic  level.  The  model  consists  of  balance  equations  for  each  extensive 
quantity  (mass,  momentum  and  energy)  that  is  transported,  constitutive  relations 
describing  the  properties  of  the  particular  phases  involved,  source  functions  of  the 
extensive  quantities,  and  initial  and  boundary  conditions  stated  at  the  macroscopic  level. 

For  descriptions  of  the  hydro-thermal  behavior  of  a multiphase  porous  medium, 
constitutive  models  for  the  constituents  are  discussed  using  macroscopic  field  variables 
such  as  permeability  and  conductivity.  For  quantitative  solution  of  the  practical 
engineering  problems  considered  in  this  study,  assumptions  made  in  the  mathematical 
model  are  discussed.  Also,  constitutive  relations  employed  within  governing  field  laws 
and  source  functions  for  field  variables  are  presented  with  mathematical  derivation  of  the 
governing  equations.  The  nature  of  boundary  conditions  relating  to  transport  phenomena 
of  the  extensive  quantities  at  the  macroscopic  level  is  briefly  discussed.  In  the  following 
chapter  when  a numerical  model  is  presented,  additional  details  of  modeling  the 
boundary  conditions  will  be  given. 

3.1  Concrete  as  a Porous  Medium 

A common  characteristic  of  concrete  as  a porous  medium  is  that  the  solid  phase  is 
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distributed  throughout  its  domain  (and  thus,  also  the  fluid  phases).  That  is,  if 
sufficiently  large  samples  of  the  concrete,  which  are  small  enough  so  as  to  represent  a 
sufficiently  close  neighborhood  around  the  point  of  sampling,  are  taken  at  different 
locations  within  a domain,  a solid  phase  is  found  in  each  of  them. 

Denoting  the  volume  of  a sample  that  satisfies  these  conditions  as  a 
representative  elementary  volume  (REV)  (shown  in  Figure  3-1)  of  the  concrete 
domain  at  the  given  point,  Bear  (1972)  and  Bear  and  Bachmat  (1998)  defined 
concrete  as  a multiphase  material  characterized  by  the  following  features 

Representative  Elementary  Volume  (REV).  A representative  elementary 
volume  can  be  determined  such  that,  regardless  of  where  we  place  it  within  a concrete 
domain,  it  will  always  contain  all  the  three  phases  (i.e.,  the  solid,  the  liquid,  and  the 
gaseous  phases  shown  in  Figure  3-1).  If  such  REV  cannot  be  found  for  a given 
domain,  the  domain  cannot  be  qualified  as  a porous  medium  domain. 

Size  of  the  REV.  The  size  of  the  REV  is  such  that  parameters  that  represent  the 
characteristics  of  the  distributions  of  the  void  space  occupied  by  the  fluid  phases  and 
of  the  solid  matrix  within  it  are  statistically  meaningful.  That  is,  the  macroscopic 
effects  of  the  microscopic  configuration  of  interphase  boundaries  and  the  actual 
variation  of  extensive  quantities  (mass,  momentum,  and  energy)  within  each  phase 
are  still  retained  in  the  form  of  coefficients,  whose  structure  and  relationship  to  the 
statistical  properties  of  the  phase  can  be  determined.  The  numerical  values  of  these 
coefficients  must  be  experimentally  measurable  in  the  laboratory,  or  in  the  field. 

Therefore,  an  REV  of  concrete  involves  (1)  a solid  phase  of  aggregates, 
hydrated  cement  and  their  bonds,  (2)  a liquid  phase  of  free  water  in  pores  and 
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chemically  bound  water  in  the  hydrated  cement,  and  (3)  a gaseous  phase  of  a mixture 
of  dry  air  and  water  vapor. 


Liquid  phase 
(liquid  water) 


Solid  phase 


dV,  REV 


Gaseous  phase 
(air-vapor  mixture) 


Figure  3-1  An  REV  occupied  by  the  three  phases 

As  stated  earlier,  the  goal  is  to  describe  each  phase  and  its  behavior  as  a 
continuous  body  that  occupies  the  entire  domain  of  the  system.  To  achieve  this 
objective,  averaging  rules  that  transform  each  of  these  phases  into  a continuum  are 
presented  in  the  following  section. 

3.2  Averaging  Principles 

Having  introduced  a representative  elementary  volume  (REV)  as  the  first  step 
necessary  for  the  process  of  integrating  from  the  microscopic  level  to  a scale  of 
continuum  mechanics,  we  now  review  the  averaging  principles  that  are  crucial  for  the 


development. 
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3.2.1  Volume  averaging  process 

First,  we  choose  a point  of  the  concrete  system  shown  in  Figure  3-1,  assuming 
that  point  to  be  the  centroid  of  averaging  volume  element  denoted  by  dV  (or  an 
REV).  Likewise,  dVs , dV"  and  dVg  will  denote  the  corresponding  subdomains  and 
averaging  volumes  of  the  solid  phase  (superscript  s ),  the  liquid  water  phase 
(superscript  w),  and  the  gaseous  phase  (superscript  g),  respectively.  Thus,  the 
averaging  volume  dV  is  the  sum  of  the  volumes  occupied  by  three  phases 


The  position  of  center  of  an  REV  in  a global  coordinate  system  is  described  by 
position  vector  x . Figure  3-2  shows  the  position  vector  denoted  by  xlocal  of  a point 
that  is  located  within  a domain  d V centered  at  x . The  volume  of  constituent  n 
within  the  REV,  d.Vn , is  obtained  by  introducing  a phase  distribution  function,  yn 


where  the  integration  refers  to  the  microscopic  volume  element  dVm  that  its  position 
is  described  by  a vector  r between  the  two  points  where  r = x'°“'  - x . 


dV  = dVs  +dVw  + dVg 


(3.1) 


for  r e d Vs  where  n ^ C, 


for  r e dVn 


(3.2) 


Thus 


<tv'(x)=  l r’(r)dVm 


(3.3) 
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z 


dVm 


Figure  3-2  Definition  of  averaging  volume  element  of  a porous  medium  (after  Lewis 

and  Schrefler,  2000) 


Similarly,  the  part  daK  of  area  da  of  the  REV,  occupied  by  the  constituent  n can  be 
written  as 


^(X)=  l/*(T)dan, 


(3-4) 


where  dam  is  the  microscopic  area  element. 

Therefore  the  volume  fraction  of  phase  n within  the  REV  becomes 


*7*  (x) 


dVn 

dV 


dV 


Lr'UJV. 


(3.5) 


which  leads  to 
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£>'=1  (16) 

K= 1 

Let  us  apply  the  volume  fraction  concept  to  an  intrinsic  microscopic  conserved 
quantity  (y/)n , which  is  averaged  over  the  entire  volume  dV  of  an  REV  at  time  t as 


^(x,t) 


r 


dV 


l r*(r,t)dvm 


(3.2) 


where  represents  the  macroscopically  averaged  value  of  the  quantity  over  the 
entire  volume  dV  of  the  REV.  Then  substituting  the  volume  fraction  r/x  (x,t)  of  a 
phase  n with  respect  to  the  entire  volume  dV  of  the  REV  into  Eqn.  (3.7),  we  have 


Tjr(x,t)=i7'(x,t)  (y/)K  (x,t)  ^ 

This  is  known  as  volume  averaging  process  (Lewis  and  Schrefler  2000),  which  will 
be  used  in  derivation  of  macroscopic  governing  equations  later. 


3.2.2  The  size  of  an  REV 

In  general,  if  an  arbitrary  elementary  volume  is  selected  as  an  REV,  then 
different  volumes  will  yield  different  averaged  values  of  each  quantity  of  interest. 
There  is  no  sense  in  asking  which  of  them  is  more  ‘correct’  unless  the  model’s 
objectives  in  a particular  case  are  clearly  stated.  Also,  in  a field  test  or  laboratory,  the 
measured  values  should  correspond  to  that  of  the  selected  arbitrary  elementary 
volume.  The  predicted  and  measured  values,  then,  can  be  the  same  within  a range  of 
error  introduced  by  the  concept  of  the  averaging  process.  However,  this  approach  has 
the  main  drawback  that  since  every  averaged  value  depends  on  the  selected  size  of  the 
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arbitrary  elementary  volume,  it  must  be  specified  by  the  size  of  the  volume  over 
which  it  was  taken.  To  eliminate  any  size  dependency  of  the  values  of  quantities,  it  is 
desirable  to  have  a universal  criterion,  which  is  based  on  measurable  characteristics 
of  porous  media. 

Using  the  criterion  for  concrete,  we  can  determine  a range  of  averaging  volumes 
within  which  the  characteristics  of  a concrete  remain  constant  (theoretically).  Then, 
computed  and  observed  values  will  be  close  enough  to  consider  an  averaged  volume 
in  that  range  as  an  REV.  Average  quantities  are  also  to  be  continuous  in  space  and 
time.  In  other  words,  the  REV  must  be  small  enough  to  be  thought  as  infinitesimal 
(i.e.  the  partial  derivatives  in  the  governing  equations  must  make  sense)  yet  large 
enough  to  give  average  quantities  without  fluctuation  depending  on  the  size  of  the 
REV. 

Figure  3-3  illustrates  the  requirements  for  the  range  of  the  REV  in  variation  of 
void  ratios  with  respect  to  averaging  volumes.  For  very  small  values  of  dV(x ) , the 
void  ratio  can  be  one  or  zero  since  x happens  to  fall  in  either  the  void  space  or  the 
solid  phase.  However,  as  the  size  of  dV(x)  increases,  large  fluctuations  in  the 
porosity  can  be  observed.  As  the  size  of  dV (x)  continues  to  increase,  we  notice  that 
the  fluctuations  gradually  decay.  When  dV(x ) = dVmm  , the  void  ratio  becomes,  more 
or  less,  constant  with  only  small  fluctuations.  The  REV  is  the  volume  within  the 
range  of  dVmm  < dV  < dVmaii  that  makes  the  void  ratio  independent  of  the  volume 
size  and  thus,  a single  valued  function  of  x only.  Then  the  void  ratio  represents  the 
porosity  of  the  porous  medium. 
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Figure  3-3  Variation  of  the  void  ratio  as  a function  of  averaging  volume  (after  Bear 

and  Bachmat  1 998) 


Considering  multiphase  mass  and  heat  transfer  in  concrete,  where  the  state  of 
each  phase  is  described  by  a set  of  variables  (permeability,  density,  and  conductivity, 
etc.),  we  undertake  a similar  approach  that  associates  with  the  geometrical 
characteristics  of  the  void  space,  e.g.  the  porosity,  with  respect  to  these  variables.  A 
range  for  REV  is  found  and  selected  for  each  variable.  Then  a common  range  of 
REV,  which  is  determined  by  a characteristic  length  k of  the  average  volume  that 
depends  on  the  specific  constituent  (or  phase)  of  the  domain,  can  be  found  for  all  of 
them.  The  characteristic  length  of  the  REV  is  required  that 
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4n*  <*<L 

where  dmm  is  the  maximum  dimension  of  a pore  or  a grain  of  concrete  and  L is  the 
characteristic  length  of  the  entire  domain  of  the  concrete.  Therefore,  for  every 
parameter  and  state  variable  of  a phase  within  a given  concrete  domain,  we  are  able  to 
define  an  average  quantity.  One  such  parameter  is  porosity  and  we  now  look  into  a 
mathematical  representation  of  porosity  by  the  volume  fractions  of  concrete. 


3,2.3  Mathematical  representation  of  the  averaged  porosity 

The  total  void  space  within  the  REV  is  defined  as  the  volume  of  the 

interconnected  void  space  dVvc  and  the  total  volume  of  the  isolated  simply  connected 
voids  dV ” that  are  occupied  by  both  the  fluid  phases 

dV*  = dVvc+dVvi  = dVw  + dVg  (39 

where  dV * represents  the  total  void  space  within  the  REV  and  dV"  and  dVg 
represent  the  spaces  occupied  by  the  liquid  water  and  the  gaseous  phases, 
respectively.  The  ratio  is  defined  as  the  total  (volumetric)  porosity  of  the  porous 
medium  at  point  x . 


n(x)  = 


dVw  + dVg 
dV 

x 


dV« 

dV 

x 


(3.10) 


From  the  point  of  view  of  fluid  flow,  only  the  volume  of  dVvc  is  referred  to  as  the 
volume  of  the  effective  void  space  subdomain.  Then  the  effective  porosity  of  the 
porous  medium  at  x is  defined  as 
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«e(x)  = 


(3.11) 


However,  throughout  this  study,  we  assume  concrete  as  a porous  medium  in  which  all 
the  void  spaces  are  interconnected.  Its  volume  is  referred  to  as  dV* , instead  of  dVvc , 
and  the  symbol  n will  be  used  to  denote  the  porosity. 

With  the  volume  fraction  v'r(x,t) , the  volume  averaging  theory  is  now  applied 
to  the  geometrical  characteristics  of  the  volume  fractions  of  the  constituents  (i.e., 
solid,  liquid,  and  gaseous  phase) 


, , (dVw  + dVg) 

dV 

w dVw 

77  =n 

( dVw  + dVg ) 

g dVg 

rjg  - n 

( dVw  + dVg ) 


= 1-77 


= n-S 


w 


= n-S 


g 


(3.12) 


with 


sw+sg=  1 


(3.13) 


where  Sw  and  Sg  are  the  degree  of  liquid  phase  (i.e.,  water)  saturation  and  the  degree 

of  gaseous  phase  (i.e.,  a mixture  of  air  and  vapor)  saturation,  respectively. 

To  summarize  the  discussion  this  far,  the  values  of  the  parameters  that  describe 
all  averaged  geometrical  characteristics  of  the  microstructure  of  the  porous  medium  at 
any  point  of  the  macroscopic  domain  are  required  to  be  single  value  functions  of  the 
location  of  that  point  and  of  time,  which  is  independent  of  the  size  of  the  averaged 
volume.  By  using  this  postulate,  particular  macroscopic  governing  equations  for 
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extensive  quantities  of  interest  are  derived  by  using  the  averaging  theory  that  has 
been  established  on  spatial  averaging  processes  in  this  section.  First,  we  will  develop 
the  generic  balance  equation  of  an  extensive  quantity  at  the  microscopic  level  in  the 
subsequent  section. 


The  formulation  of  a generic  balance  equation  for  an  extensive  quantity  (e.g., 
energy)  in  the  neighborhood  of  an  arbitrary  point  within  continua  for  single  materials 
is  presented.  Consider  a point,  which  is  located  within  a fixed  finite  domain,  U, 
referred  to  as  a control  volume  and  bounded  by  a closed  surface  dU  as  shown  in 
Figure  3-4. 


Figure  3-4  A control  volume  U bounded  by  a closed  surface  dU  and  outward  unit 

normal  vector  n 

For  a generic  conserved  variable  y/  defined  per  unit  mass  and  describing  the 
microscopic  situation  of  the  body,  the  general  form  of  a balance  law  over  U of  the 
body  may  be  expressed  as 


3.3  Balance  Laws  for  Continua 


n 


cU 


l 


rate  change  of 
y/  within  U 


- {the  flux  of  y/  across  dU } - 


(3.14) 
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In  mathematical  form,  the  balance  components  of  Eqn.  (3.14)  can  be  expressed  as 
(a)  the  rate  change  of  the  amount  of  y/  within  U, 


H[P¥da  <315> 

dt  >' 

(b)  the  flux,  i , into  U through  dU  , 

Cf  i.n  dS  <316> 

Jeu 

and  (c)  the  external  supply,  g , of  y/ , 


[pgdU 


(3.17) 


where  p is  the  mass  density  of  y/  (i.e.,  the  amount  of  y/  per  unit  volume  of  the  body) 
and  n is  outward  unit  vector  normal  to  dU  . 

Substituting  the  mathematical  expressions  above  into  Eqn.  (3.14),  we  now  have 
the  general  form  of  a balance  equation  of  y/  within  U 


jtlPWdn-^j.nds-[pgda=o  (318) 

For  example,  if  we  set  y/  = 1 , i = 0 , and  g - 0 , then  Eqn.  (3.18)  becomes  the 
mass  balance  equation  of  a simple  continuum 

H[pdn= o <319> 

dt  i l 

which  states  that  the  total  mass  inside  a material  volume  remains  constant  for  a 


continuum.  For  the  momentum  balance,  the  momentum  per  unit  mass  or  velocity 
v(x,t)  is  taken  for  y/  while  we  set  i = a , and  g = b . Then,  we  have 
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1 pv  dQ  - ^ cr»n  dS  - | pb  dQ  = 0 


(3.20) 


where  the  momentum  flux  is  represented  by  cr  (stress)  and  the  body  force  per  unit 
mass  acting  on  the  body  is  denoted  by  b (the  gravitational  acceleration). 

If  y/  is  used  to  represent  the  internal  energy  E plus  the  kinetic  energy  ^ vv 


per  unit  mass,  y/  = E + ^-vv , the  energy  balance  in  the  body  can  be  expressed 


d 

dt 


J p{E  + — v»v)  dQ.  + cj^  (q  - cr*v)»n  dS  - J (h  + b»\)  dQ.  = 0 


(3.21) 


where  the  energy  flux  has  two  parts  i = -q  + <r»v , that  is,  -q  represents  the  heat  flux 
and  <r»v  accounts  for  the  work  done  by  stress.  The  external-energy  supply  term,  g , 
in  Eqn.  (3.18)  is  replaced  by  the  rate  of  external  heating  (e.g.,  radiation),  h,  and  the 
work  done  by  body  forces,  b»\ . 


3.3.1  Governing  differential  equations  for  microscopic  conserved  quantities 

Recall  the  following  two  theorems  (1)  Gauss  theorem,  0 f*n  dS  = [ V*f  dQ 


for  a differentiable  vector-valued  function  f(x,t)  on  a region  containing  U and  8U , 


and  (2)  the  Reynolds  transport  theorem,  — f f dQ  = [ 

dt  J •* 


[ f dQ  = [ 

^ + V-(v  f) 

JU  JU 

L dt 

dQ , where 


v(x,/)  is  the  velocity  of  the  material. 

Applying  the  Reynolds  transport  theorem  to  the  first  term  in  Eqn.  (3.18)  and  the 
Gauss  theorem  to  the  second  term,  we  have 
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■ d(PV) 
,J  dt 


+ wV(py)  + py/W»\ -V»i- pg 


d Q = 0 


(3.22) 


However,  if  f(x,/)  is  a continuous  real-valued  function  defined  on  a region 
containing  U and  dU  and  j f(x,t)  dQ.  = 0 for  every  subregion  of  U,  then  f(x,/) 

vanishes  identically  on  U,  i.e.,  the  fundamental  lemma  of  calculus  of  variations. 
Therefore,  the  terms  in  the  bracket  in  Eqn.  (3.22)  can  be  written 


+ yS7(py)  + py  V*v  - V»i  - pg  - 0 
dt 


(3.23) 


for  any  arbitrary  subregion  of  U.  By  defining  the  material  derivative  operator 
Dt  Dt  = 8/  dt  + v*V , which  measures  the  time  rate  of  change  at  a specific  material 
point,  the  local  balance  equation  can  be  expressed 


D(pip)  , w „ A (3.24) 

— CJ—  + py/  V»v-  V»i-yC>e  = 0 v ’ 

Dt 

Therefore,  if  we  set  y - 1 , i = 0 , and  g = 0 , the  local  mass  balance  law 
defined  by  Eqn.  (3.19)  is  rewritten  as 


RP.  + pV.v  = o 

Dt 

For  the  momentum  the  local  balance  law  is  rewritten 


(3.25) 


D(pv)  + pVv.v-  V»<j  - pb  = 0 
Dt 


(3.26) 


Since  D(pv)t  Dt  = pDv  / Dt + vDpt  Dt , 
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(3.27) 


However,  by  virtue  of  the  mass  balance,  the  term  in  parentheses  in  the  last  equation 
vanishes.  Thus  we  have 

Dv  ^ (3.2 

p V»cr-  pb  = 0 v 

Dt 

which  is  known  as  the  linear  momentum  balance  law  or  Cauchy’s  first  law  of 
continuum  mechanics. 

For  the  local  energy  balance,  Eqn.  (3.21)  is  rewritten  in  the  form  of 


If  we  expand  the  material  derivative  and  use  differential  identity  such  as 


V»(cr»v)  = v»(V»cr)  + cr  : Vv 


(3.30) 


then  we  have 


y Dt  A 2 ) 


(3.31) 


By  eliminating  the  terms  that  vanish  in  accordance  with  the  mass  balance  in  the  first 
parenthesis  and  by  forming  the  dot  product  of  v with  both  sides  of  the  momentum 
balance  in  the  last  parenthesis 
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(3.32) 


and  the  energy  balance  law  can  be  expressed  in  the  form  of 


p + V»q-cr  : Vv-  ph  = 0 


(3.33) 


Dt 


The  terms  in  the  equation  (3.33)  have  the  following  physical  interpretations 
pDE  / Dt  is  the  rate  of  change  of  internal  energy,  V*q  is  the  rate  of  heat  flow, 

- <j  : Vv  represents  heating  caused  by  compression  and  dissipative  momentum 
transfer  (e.g.,  plastic  deformation),  and  - ph  is  the  rate  of  heating  by  external 
supplies. 

3.3.2  Governing  differential  equations  for  multi-phase  continua 

The  balance  laws  for  a multi-phase  continuum  are  analogous  to  those  of  a 
simple  continuum,  except  that  global  balances  apply  to  the  multi-phase  mixture  as  a 
whole  via  an  averaging  process  used  for  a frame  work  to  construct  mathematical 
formulations  of  governing  field  equations  at  the  macroscopic  level. 

To  start  with,  the  general  balance  law  for  a generic  extensive  quantity  defined  in 
Eqn.  (3.24)  is  rewritten  for  a multiphase  mixture,  which  is  a collection  of  overlapping 
continua  called  phases  (or  constituents) 


where  each  constituent  n occupies  a fraction  of  the  overall  mixture  volume  defined 


y Dn  (77* p*y/*) 

rw 


1 .ft  v7  , Y7  *-ji  r\ 

+ Tj  p yp  V»v  -V'7  1 -77  p g = 0 


(3.34) 


by  a volume  fraction  7?r(x,t) , and  has  the  intrinsic  mass  density  p * with  its  own 
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velocity  v*(x,t).  Introducing  the  material  time  derivative  of  any  differential  function 
f*(x,t)  of  the  k phase  as 


DKfK 

Dt 


dfn 

— — h grad  f*  •v* 
dt 


(3.35) 


we  can  define  the  time  derivative  of  /*(x,/)  taken  with  respect  to  another  moving 
phase  a 


where 


d°  r 

Dt 


dr 

— — + grad /*  •\an 
dt 


(3.36) 


yan  = v“  -v* 


(3.37) 


While  the  equation  (3.34)  provides  balance  laws  for  the  overall  mixture,  it  can 
be  divided  for  each  constituent’s  balance  laws 


D\ 7?W) 

Dt 


i y”7  , .ft  V7  ~*ft\ft  ~.ft  ~ft  CLft  „ft 

+ rj  p \y  V*v  —\»rj  l —rj  p g —e 


(3.38) 


The  term  en  on  the  right  hand  side  accounts  for  all  exchanges  of  the  conserved 
quantity  y/K  into  constituent  n from  other  constituents.  It  should  be  noted  that,  in 
order  to  conserve  the  balance  of  the  overall  multi-phase  mixture,  we  must  have 


Y.e*=  0 


K 


(3.39) 
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3.4  Macroscopic  Mass  Balance  Law  for  a Multiphase  Mixture 

From  Eqn.  (3.38),  we  will  derive  the  macroscopic  mass  balance  law  for  each 
individual  constituent  denoted  by  n (i.e.,  the  solid,  the  liquid,  and  the  gaseous  phase) 
in  a multiphase  mixture. 

3.4.1  The  mass  balance  law  for  the  solid  phase 

As  with  the  mass  balance  for  the  solid  phase  (superscript  s)  as  a single 

continuum,  we  set  = 1 , \n  = 0 , and  g”  = 0 in  Eqn.  (3.38) 


°nP  +7 fps  div(vs)  = 0 
Dt 


(3.40) 


However,  since  77s  =1-/7  from  Eqn.  (3.12),  the  macroscopic  mass  balance  law  for 
the  solid  phase  is  rewritten 


Z)S[(1-«)PS1 

— 4- + (1  -n)  ps  div  (■ vs ) = 0 


or,  since  D/ Dt  = d/ dt  + v«V , 


(3.41) 


d\(\-n)ps]  r , (3.42 

— - — - + div[(l-«)/?s  vsJ  = 0 

where  (1  -n)ps  is  the  macroscopic  volume-averaged  mass  density  of  the  solid  phase 
and  vs  is  its  velocity.  It  is  assumed  that  there  is  no  mass  exchange  (e.g.,  chemical 


interaction)  between  the  solid  phase  and  the  other  phases,  i.e.,  es  = 0 . 
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By  using  the  vector  identity, 

div  [(1  - n)ps  vs  J = (1  - n ) ps  div  vs  + grad  [(1  - «) ps  ]»vs 
and  assuming  grad[(l  -n)/?5]  = 0 , Eqn.  (3.42)  becomes 


d\(\-ri)p  1 

L J+(l  n)ps  div  vs  = 0 

dt 

(3.44) 

Expanding  and  then  dividing  both  sides  by  p% , we  have  another  form  of  the 
macroscopic  mass  balance  of  the  solid  phase 

(l-«)5/7s  dn  s n 

+ (l-«)divv  =0 

ps  dtdt 

(3.45) 

or 

(1  -n)Dsps  D'n  n ...  s _ 

— + (l-«)  divv  =0 

ps  Dt  Dt 

(3.46) 

3.4.2  The  mass  balance  law  for  the  fluid  phases 

Again,  as  with  the  mass  balance  for  a fluid  phase  n 

(where  n ^ s')  as  a single 

continuum,  we  set  = 1 , \K  = 0 , and  g*  = 0 in  Eqn.  (3.38) 


DK  (rf  p") 
Dt 


+ rf  pK div  v"  = ±rf p*e" 


(3.47) 


where  rf  p"  is  the  macroscopic  volume-averaged  mass  density  of  the  fluid  phase, 
\n  is  its  velocity,  and  rf pKen  is  the  rate  of  mass  exchange  between  the  fluid  phases 


(e.g.,  -rf p”en  is  the  quantity  of  water  lost  through  evaporation  per  unit  time  per  unit 
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volume  and  +ti*  p^e*  is  the  quantity  of  vapor  gained  through  condensation  per  unit 
time  per  unit  volume). 

By  replacing  the  volume  fraction  r/n  by  nSn  (see  Eqn.  (3.12))  and  the  rate  of 
mass  exchange  term  rf  pnen  by  m , 


D'inS'P*) 

Dt 


+ nSr[pK div  v*  =±m 


(3.48) 


where  S K is  the  degree  of  saturation  of  the  fluid  phases. 

Therefore,  the  macroscopic  mass  balance  of  the  liquid  phase  (i.e.  water)  is 
defined 


D’inS'P*) 

Dt 


+ nSwpw  divv*  =-m 


(3.49) 


where  -m  is  the  quantity  of  water  lost  through  evaporation  per  unit  time  per  unit 
volume. 


3.4.4  Species  mass  balance  in  the  gaseous  phase 

Recall  that  a collection  of  overlapping  continua  called  phases  was  defined 
earlier  as  a multi-phase  mixture.  The  constituents  (or  phases)  can  be  segregated  at  the 
microscopic  level,  where  the  phases  behave  as  continua  but  their  small-scale  motions 
are  inaccessible  to  field  measurement.  Such  mixtures  are  named  multi-phase 
mixtures. 

However,  the  gaseous  phase  considered  in  this  study  is  assumed  to  be  a perfect 
mixture  of  two  ideal  gas  species,  i.e.  dry  air  and  vapor.  That  is,  the  single  gas  phase 
can  be  decomposed  into  two  different  species  where  the  segregation  occurs  at  a 
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molecular  scale  but  not  at  the  microscopic  scale.  It  is  called  a multi-species  mixture, 
the  species  (or  components),  dry  air  and  vapor,  of  the  gaseous  phase  (i.e.,  moist  air), 
are  considered  to  be  miscible  and  they  occupy  the  same  volume  fraction.  That  is, 

p*  = pga  + pgw  ('3‘5< 

where  the  mass  density  of  the  gas  phase  is  the  sum  of  the  mass  densities  of  dry  air 
(superscript  ga ) and  vapor  (superscript  gw).  Introducing  the  mass  fraction  of 
component  a = ga,  gw , 


CO  — p / p 


(3.51) 


By  using  the  mass  fraction,  the  velocity  of  the  gas  phase  vg  can  be  expressed 


Vs  = — (pgavga  + p^v*")  = coga\ga  + CO ^v^  (3-52) 

pg 

This  is  called  the  barycentric  velocity,  which  is  the  mass-weighted  mean  velocity. 
Furthermore,  this  implies  that  a difference  between  the  species’  velocities  and  the 
mean  velocity  of  the  gaseous  phase  exists.  We  define  the  velocity  difference  as  the 
diffusion  velocity  of  component  a with  respect  for  the  mean  gas  flow,  u" 


From  the  macroscopic  mass  balance  of  the  fluid  phases  of  Eqn.  (3.48),  we 
introduce  the  diffusion  velocity  into  the  macroscopic  mass  balance  law  for  each 
component  of  the  gaseous  phase,  first,  for  dry  air, 


(3.53) 
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d("Sgpga ) 

dt 


-+ 


div  \jiSgpga  ugo  J + div  [nSgpga  vg  J = 0 


(3.54) 


and  next  for  vapor, 


d{nSgPn 

dt 


■dxv\jiSgpgw  u*M'J  + div[^«5,g/7*w  vg  J =m 


(3.55) 


where  the  terms  nSgpa  u"  accounts  for  the  diffusive  mass  flux  of  component  a 
By  introducing  the  diffusive  mass  flux  of  each  component  a as 


}a=nSgpaua 


(3.56) 


the  macroscopic  mass  balance  law  for  dry  air  of  the  gaseous  phase  is  then  rewritten 


d(nSgpga) 

8t 


+ div  jga  + div  ^nSgpga  vg  J = 0 


and  thus  the  macroscopic  mass  balance  law  for  vapor  becomes 


(3.57) 


d{nSgPn 

dt 


+ divj^  + div[«5,gpgw  vg]  = m 


(3.58) 


where  the  first  term  on  the  left  side  accounts  for  the  accumulation  of  component  a 
known  as  the  mass  concentration  of  a in  the  gaseous  phase,  i.e.  moist  air,  the  second 
term  models  the  diffusion  of  a with  respect  to  the  mean  macroscopic  flow  of  moist 
air  in  the  porous  medium,  the  third  term  models  the  advection  of  a with  the  mean 
macroscopic  flow  velocity,  and  the  right  hand  side  of  the  equation  accounts  for  the 
exchange  of  a to  or  from  the  other  molecular  species  in  the  phase  or  the  other  phase, 
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i.e.  the  liquid  phase  in  the  porous  medium.  The  species  mass  balance  law  is  later  used 
to  obtain  general  field  equations,  which  will  be  numerically  solved  in  Chapter  4. 

3.5  Macroscopic  Linear  Momentum  Balance  Law  for  a Multiphase  Mixture 

We  now  look  into  applications  of  mixture  theory  to  the  equations  governing 
fluid  flows  through  porous  media.  We  will  start  from  the  generic  balance  equation  of 
a single  constituent,  Eqn.  (3.38)  to  derive  the  momentum  balance  for  a single  fluid 
phase  flowing  through  a concrete  matrix.  With  the  concept  of  the  volume  fraction,  the 
macroscopic  momentum  balance  equations  of  the  fluid  phases  will  be  extended  to 
Darcy’s  law,  the  fundamental  field  equation  for  porous-medium  flow. 

3.5.1  The  linear  momentum  balance  law  for  the  solid  phase 

If  we  set  \f/"  - v* , i*  = ct*  , gK  =b* , and  eK  = r*  (i.e.,  the  rate  of  momentum 
exchange  into  the  phase  tc)  in  Eqn.  (3.38),  then  we  have 


Dt 


+ TjK  PKV*  V.V*  - v.  rfoK  - T]K  pKbn  = rn 


(3.59) 


Eliminating  terms  on  the  left  side  of  this  equation  in  a manner  similar  to  that  used  in 
derivation  of  the  momentum  balance  law  of  a simple  continuum  (see  also  the 
derivation  of  Eqn.  (3.28)),  we  have  the  macroscopic  linear  momentum  balance  law 
written  as 


DK\n 

rf  pK  ^ — V*  rf  a*  - rf  pfb’  = r* 


(3.60) 


or 
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dt 


(3.61) 


Let  us  neglect  the  acceleration  terms  as  a first  approximation  for  slow  flow  of  a 
macroscopically  averaged  fluid  phase  through  a porous  medium.  Therefore,  we 
simplify  Eqn.  (3.60)  for  the  macroscopic  momentum  balance  equation 


where  bn  is  replaced  by  rf  g , which  accounts  for  the  averaged  gravitational  effects. 

Now  we  turn  to  the  momentum  balance  law  for  the  solid  phase.  Neglecting  the 
velocity  of  the  solid  phase,  vs  = 0 , we  further  simplify  Eqn.  (3.60)  for  the  linear 
momentum  balance  equation  of  the  solid  phase 


where  crs  represents  the  stress  tensor  in  the  solid  phase,  g represents  the  gravitational 
acceleration  , and  rs  represents  the  momentum  exchange  to  the  other  fluid  phases. 

3.5.2  The  linear  momentum  balance  law  for  the  fluid  phases 

For  the  fluid  phases,  we  begin  with  two  assumptions  (1)  The  fluid’s  inertial 
effect  is  assumed  to  be  negligible  on  its  momentum  balance 


(3.62) 


(3.63) 


where  n = w,  g 


(3.64) 


Dt 
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and  (2)  the  gravity  is  the  only  body  force  in  the  form  of 

Li  , (3.65 

b = g where  n = w,g 

Therefore,  we  have  the  linear  momentum  balance  equation  for  the  fluid  phases 
written  as 

V'tfcr"  +rj7tp’Tg  = —r*  where  7t  = w,g 
Since  vs  = 0 , the  momentum  exchange  rK  of  the  fluid  phases  to  the  solid  phase 
can  be  expressed  by  Stokes  ’ drag  (Allen  et  al.  1988) 


where  the  constant  yK  is  called  the  fluid  mobility.  Then  the  linear  momentum 
balance  equation  of  the  fluid  phases  becomes 


In  addition,  the  counterpart  of  the  momentum  exchange  of  the  solid  phase  to  the 
fluid  phases  in  Eqn.  (3.63)  becomes 


(3.67) 


(3.68) 


(3.69) 


in  order  to  meet  the  requirement  in  Eqn.  (3.39). 


3.5.3  The  linear  momentum  balance  law  for  the  multiphase  medium 


Through  summing  all  the  momentum  balance  equations  written  for  the  solid  and 
fluid  phases  and  by  assuming  continuity  of  stress  field  across  the  solid-fluid  interface. 
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the  linear  momentum  balance  equation  for  the  multiphase  medium  as  a whole  is 
obtained  as 


V.  (77  V*  + rfcj* ) + (7>v  + rf  p*  )g  - if  pK  ( v*  .grad  v* ) = 0 


(3.70) 


where  n = w,  g . 

Further  simplification  can  be  made  on  the  last  equation  by  neglecting  the  last 
term  on  the  left  hand  side  that  depends  on  the  gradient  of  the  fluid  velocity,  by 
defining  the  total  stress  acting  on  a unit  area  of  the  multiphase  medium 


erto,a/  = 77  Vs  + rjwaw  + rfo*  = (1  - n)as  + n(Swaw  + Sng ) (3'71) 


and  by  introducing  the  averaged  mass  density  of  the  multiphase  medium 


pmg  = 7 fp5  +7 fpw  + r1gpg  = (1  - n)ps  + n(Swpw  +S  ps)  (3j2) 


Therefore  the  linear  momentum  balance  equation  of  the  multiphase  medium  becomes 


S7^(J,0,a,  + pmgg  = 0 


(3.73) 


where  cr'°,al  is  the  total  stress  acting  on  a unit  area  of  the  multiphase  medium  and 


pmg  is  the  volume-averaged  mass  density  of  the  multiphase  medium. 


3.6  Macroscopic  Energy  Balance  Law  for  a Multiphase  Mixture 

Again,  we  shall  begin  with  the  generic  macroscopic  balance  equation  of  Eqn. 
(3.38)  by  taking  the  following  components  into  account  in  the  generic  equation.  First, 

if  we  set  \pK  - E"  + — vK  •v’1 , then  the  generic  extensive  quantity  if/"  is  considered  as 
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the  internal  energy  per  unit  mass  E n plus  the  kinetic  energy  per  unit  mass  — v^v* . 

The  flux  of  energy  \n  has  also  two  parts  we  set  i"  = t]k{gk  »v*  -q*)  where  cr'^v* 
accounts  for  the  work  done  by  stress  and  - q*  represents  the  heat  flux.  For  external 
heating  supplies,  we  have  the  rate  of  external  heating  (e.g.,  attributable  to  radiation) 
hn  and  work-done  by  body  forces  denoted  by  b"»\n . Finally,  the  exchange  of  energy 
between  phases  due  to  phase  change  is  substituted  for  en . To  summarize,  we  have  the 
following  expressions 


y/*  =E"+  -v'.v* 

2 

i'  = 17,(0^  •v'-  q')  (3-74) 

gK  = hn+b’c*\n 
e * = pnRK 

Therefore,  substituting  these  components  into  Eqn.  (3.38),  we  have 


Dn  1 f p*{EK  +-vK»v*) 


Dt 


+ rj*  p*(E*  + -Iv*»v*)V»v’r 


(3.75) 


- V-  -v*  - q' ) - rfp*  (i h?  + b * .v' ) = pnRn 


By  a similar  approach  to  simplify  the  formulation  of  Eqn.  (3.33) — the  generic 
local  energy  balance  equation— Eqn.  (3.75)  can  be  written  in  a simplified  form  as 


DnFn 

77  p — — YTj  V»q  —77  cr  :Vv  —77  p (n  + g»v  ) = tj  p K 


(3.76) 


or 
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7'  | •(/>'£')  + n^-(p’E- ) + 7”  V-q'  (3.77) 

- 77  V"  : Vv"  - 7 V"  (/?*  + g»\* ) = rf  pnRK 

where  E * is  the  intrinsic  specific  energy  of  a control  volume  of  a phase  n , cf  is  the 

internal  heat  flux  in  the  phase,  h”  is  the  heat  source,  and  R*  is  the  exchange  of 
internal  energy  (e.g.  latent  heat)  between  phases  due  to  phase  change.  A physical 
interpretation  of  the  last  equation  can  be  given  in  that  the  first  term  on  the  left 
represents  the  rate  of  change  of  internal  energy,  V*q'T  represents  the  rate  of  heat  flow, 
- a*  : Vv*  accounts  for  heating  by  mechanical  deformation  (e.g.,  plastic  deformation 
or  viscous  friction),  -rf  p”  (hK  +g»v'r)  accounts  for  external  heating  supplies  plus 
external  work  done  by  body  forces,  and  rj*  p*R*  accounts  for  energy  exchange  of  a 
phase  through  phase  change.  This  is  the  macroscopic  energy  balance  equation  for  a 
constituent  of  a multiphase  continuum. 

Therefore,  the  energy  balance  of  the  total  system  can  be  written  as 


z 


7 p 


D*E* 

Dt 


+ 7*V<  -7V'  : Vv*  -ri*ff(h*  + g^n)-pK pKR 


= 0 


(3.78) 


which  physically  means  that  the  total  balance  of  energy  exchange  among  all  the 
phases  is  zero.  After  introducing  constitutive  relationships,  we  will  come  back  to 
refine  the  macroscopic  energy  balance  equation. 


3.7  Constitutive  Equations 

While  the  formulation  of  governing  equations  sketched  so  far  provides  a basis 


for  the  analysis  of  continuum  phenomena  of  mass  and  heat  transfer  in  a porous 
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medium,  they  do  not  contain  information  regarding  the  properties  (e.g.  mechanical 
properties  of  a particular  continuum),  which  are  necessary  to  form  a determined 
system. 

The  additional  information  is  contained  in  what  are  known  as  constitutive 
equations.  The  form  of  relationships  between  the  fluxes  and  driving  forces  in  the 
mass-,  momentum-  and  energy-balance  equations  is  usually  taken  by  constitutive 
equations.  The  linear  relationship  between  stress  and  strain  in  an  elastic  solid  is  an 
example  of  a constitutive  assumption.  Other  examples  include  that  the  thermal  flux  is 
a linear  function  of  the  temperature  gradient  and  the  mass  diffusion  flux  of  a species 
in  a binary  mixture  of  an  ideal  gas  is  linearly  proportional  to  the  concentration 
gradient  of  the  species. 

Despite  the  fact  that  these  constitutive  equations  are  often  based  on  assumptions 
used  for  defining  the  behavior  of  ideal  continua  to  construct  mathematical  models, 
they  are  arrived  at  by  physical  observation  and  proved  by  experimental  evidence. 
Therefore,  our  aim  in  this  section  is  oriented  to  develop  a conceptual  model  in  which 
simplification  is  carried  out  to  the  point  where  the  mathematical  model  is  amenable, 
yet  not  so  crude  as  to  miss  the  features  of  the  phenomena  it  is  intended  to  describe. 

For  quantitative  solution  applicable  to  practical  engineering  problems,  we 
choose  constitutive  equations  that  are  based  upon  quantities  measurable  in  the 
laboratory,  and  that  have  been  extensively  validated  both  with  reference  to  exact 
solutions  and  to  experiments.  We  will  begin  with  constitutive  equations  that  describe 
the  material  properties  of  the  fluid  phases  and  only  briefly  deal  with  the  solid  phase. 
Chapter  5 will  deal  with  this  subject  for  the  solid  phase  in  detail,  where  a stress 
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analysis  of  a multiphase  continuum  under  thermal  loading  conditions  is  the  primary 
focus. 


3.7.1  Flow  equations  in  porous  media 

Having  established  the  equation  governing  fluid  flows  through  porous  media, 
namely,  the  macroscopic  linear  momentum  balance  equation,  we  now  turn  to  a 
specific  application  to  the  momentum  balance  for  a single  fluid  phase  flowing 
through  concrete.  First,  we  start  by  deriving  Darcy’s  law,  the  fundamental  field 
equation  for  porous-media  flow.  Then  we  will  extend  Darcy’s  law  to  multi-fluid 
flows  and  present  constitutive  field  equations  applicable  to  fluid  flows  through 
concrete. 

It  is  assumed  that  each  fluid  phase  considered  in  this  study  is  Newtonian  (i.e., 
momentum  transfer  via  shear  stresses  is  negligible  compared  with  momentum 
exchange  to  the  solid  matrix).  Thus,  we  can  write 


a*  - -pn I for  n * s \ ) 

where  px  presents  the  hydrodynamic  pressure  of  the  fluid  phase  and  I stands  for  the 
identity  tensor. 

Second,  we  assume  that  gravity  is  the  only  body  force  and  has  the  form  of 


bK  = gVz  for  7r  * s 


(3.80) 


where  g denotes  the  gravitational  acceleration  and  z measures  depth  below  a datum 


level.  Applying  these  assumptions  to  Eqn.  (3.68),  we  find  that  the  fluid’s  linear 


momentum  balance  can  be  written  as 
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T1  KyK -~yn  {rfNpn-r]K  pK  gNz)  where  n = w,g 

Introducing  the  macroscopic  fluid  mobility  as 


(3.81) 


rfr 


n 


(3.82) 


where  represents  the  dynamic  viscosity  of  the  fluid  and  has  units 
[(mass)(length)~l  {time)~'  ] and  K is  the  intrinsic  permeability  [(length)2]  of  the  solid 
matrix,  we  can  rewrite  the  momentum  balance  equation  for  the  fluid  as 


i’''*=~<yp,-p’gvz)  <383) 

This  is  Darcy ’s  law,  the  field  equation  for  the  macroscopic  averaged  velocity  of  the 
fluid  of  a single  fluid  flow  in  a porous  medium. 

Despite  the  fact  that  the  permeability  K in  Eqn.  (3.83)  is  supposed  to  be  an 
intrinsic  material  property  of  the  solid  matrix,  fluid  flows  in  actual  field  experiments 
show  a remarkable  discrepancy  between  the  flow  of  air  and  the  flow  of  liquid  used  as 
a means  for  determining  permeability,  which  may  be  seen  as  an  error  in  application  of 
laminar  flow  theory  on  which  Darcy’s  law  is  based. 

3.7.2  Deviation  from  Darcy’  law 

Deviation  from  Darcy’s  law  has  been  observed  in  gas  flow  at  low  pressure  the 
flow  is  faster  than  would  be  predicted  by  Darcy’s  law.  Darcy’s  equation  was 
developed  as  a linear  relationship  between  velocity  of  the  fluid  and  the  gradient  of 
pressure,  which  is  an  approximation  to  represent  the  laminar  flow  of  Newtonian 
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liquids  through  porous  media.  However,  on  the  macroscopic  scale  of  observation, 
Darcy’s  law  does  not  accurately  predict  gas  flow  rate  through  a complex  network  of 
porous  media  such  as  concrete  where  the  pore  size  generally  being  considered  is  very 
small. 

Contrary  to  laminar  flow  theory,  on  which  Darcy’s  law  is  based  and  where  we 
assume  zero  velocity  in  the  fluid  at  the  solid-fluid  interface,  the  layer  of  gas  next  to 
the  surface  is  in  motion  along  the  solid  surface.  That  is,  if  the  solid  has  a zero 
velocity,  the  velocity  of  the  gas  layer  in  the  immediate  vicinity  of  the  wall  still  has  a 
finite  value,  not  necessarily  zero.  This  contributes  an  additional  flux  to  the  gas 
flowing  through  a porous  medium  that  is  larger  than  what  is  expected  from 
Poiseuille’s  law  (Klinkenberg  1941).  This  phenomenon  is  called  the  slip  flow  effect  or 
Klinkenberg ’s  effect. 

Klinkenberg  (1941)  demonstrated  that  the  slip  flow  phenomenon  that  arises 
during  gas  flow  (using  a glass  capillary  tube  as  a model)  was  found  to  be  quite 
substantial  for  low  permeability  materials.  He  derived  the  following  expression  for 
the  relationship  between  permeability  of  a porous  medium  to  gas  and  to  liquid 


k*=Ke 


1 + 


4 cl 
r ) 


= K( 


1 + — 


\ 

P , 


(3.84) 


where  k 8 is  the  permeability  to  the  gas,  Kc  is  the  permeability  to  liquid,  X is  the 

mean  free  path  of  the  gas  molecules  under  the  mean  pressure  p at  which  k8  is 
determined,  c is  a proportionality  factor,  r is  the  radius  of  the  capillary  tube  in 
Klinkenberg’ s model,  and  b is  a constant  for  a gas-solid  system  that  depends  on  the 
mean  free  path  of  the  gas  and  the  size  of  pore  channels  in  the  porous  medium.  Thus,  a 
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physical  interpretation  of  Klinkenberg’s  effect  can  be  that  whenever  the  mean  free 
path  of  the  gas  molecules  approaches  the  dimensions  of  the  flow  conduit,  the 
individual  gas  molecules  are  in  motion  at  the  solid-gas  interface  and  contribute  an 
additional  flux. 

By  applying  Eqn.  (3.84)  to  Eqn.  (3.83),  a generalized  Darcy’s  flow  equation  for 
gas  flow  in  porous  media  results  in 

t1sxs=_^l[  i + A 
/'H  Pgj 

3.7.3  Liquid  permeability 

In  contrast,  liquids  usually  obey  no-slip  boundary  conditions  at  the  liquid-solid 
interface.  However,  in  the  case  of  water  permeability  testing,  due  to  the  low  flow  rate, 
steady  state  flow  can  usually  be  achieved  only  after  long  periods  of  time  under  high 
pressures  that  may  cause  changes  in  the  solid  skeleton.  Besides,  in  case  of 
cementitious  materials  such  as  concrete,  these  changes  include  “continuing  hydration 
and  autogeneous  healing,  silting  of  fines,  dissolution  and  transport  of  Ca(OH)2  from 
cement  gel,  and  crystallization  in  larger  pores”  (Dhir  et  al.  1989).  This  fact  has  been 
explained  by  swelling  and  hydration  of  cement  hydrates  in  contact  with  water,  which 
influences  permeability  measured  using  liquid  water  (liquid)  flow  measurements. 

Therefore,  the  permeability  of  cementitious  materials  to  liquid  water  should  not 
be  thought  as  being  the  sole  intrinsic  parameter  of  the  solid  matrix.  To  account  for 
this  phenomenon,  a scale  factor  is  introduced  in  relation  with  the  intrinsic  gas 


(yp„-pggWz) 


(3.85) 


permeability 
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Ke=CKg  (3‘86) 

where  Kf  is  a measured  permeability  to  liquid  water,  C,  is  a factor  that  accounts  for 
changes  of  permeability  of  the  solid  matrix  to  liquid  due  to  a physical/chemical 
interaction  of  the  solid  matrix  in  contact  with  water,  and  Kg  is  the  intrinsic 

permeability  of  a cementitious,  porous  media  to  gas. 

Therefore  we  replace  Klinkenberg’s  expression  of  Eqn.  (3.84)  by  a relationship 
between  experimentally  measured  permeability  and  the  intrinsic  permeability  of  the 
cementitious,  porous  media  to  gas 


kS=Kg 


f 

V 


(3.87) 


where  Kg  is  the  intrinsic  permeability  of  a porous  medium  to  gas,  k8  is  the  gas 

permeability  measured  in  experiments,  b is  a constant  which  quantifies  the  slip  flow 
effect  (called  Klinkenberg’s  constant),  and  pg  is  the  magnitude  of  hydrodynamic 

mean  pressure  in  the  gas  flow. 

It  is  desirable  to  express  the  permeability  as  an  intrinsic  permeability  (also 
called  absolute  permeability)  for  use  in  numerical  simulations.  Neither  gas 
permeability  nor  liquid  permeability  can  be  considered  to  be  the  single  absolute 
measure  of  the  intrinsic  permeability  of  the  solid  skeleton.  This  implies  that  the 
macroscopic  permeability  of  a material  is  a function  of  both  the  solid  matrix  and  the 
particular  fluid  used  to  experimentally  measure  the  permeability. 
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For  a complete  description,  fluid  flows  in  concrete  will  be  characterized  by  two 
permeability  parameters  intrinsic  gas  permeability  Kg  , and  liquid  permeability  Kt , 

as  well  as  the  Klinkenberg’s  constant  b which  characterizes  gas  slip  flow  effects. 
These  three  parameters  taken  in  conjunction  can  be  used  to  fully  characterize  the 
single-phase  flow  of  either  liquid  or  gas  through  concrete.  The  quantification  of  these 
parameters  in  concrete  will  be  presented  in  Chapter  4 when  numerical  modeling  of 
gas  and  liquid  flow  through  heated  concrete  is  discussed. 

3.7.4  Multi-fluid  flow  equation 

While  the  discussion  presented  so  far  is  useful  for  modeling  in  single-fluid 
flows  in  porous  media,  the  flows  that  occur  in  a partially  saturated  porous  medium 
exposed  to  nonisothermal  conditions  (e.g.  high  temperatures)  involve  two  distinct 
fluid  phases  (e.g.  liquid  water  and  gaseous  steam)  flowing  simultaneously.  In  this 
subsection,  we  illustrate  the  extension  of  the  single-fluid  flow  to  multi-fluid  flows  by 
considering  a three-phase  porous  medium  whose  constituents  are  the  solid  ( s ),  the 
liquid  water  (w),  and  the  mixture  of  air  and  vapor  (g). 

Remember  that  a porous  medium  exhibits  different  permeabilities  to  different 
fluids.  This  dependence  implies  that  the  presence  of  the  liquid  phase  can  reduce  the 
effective  permeability  of  the  porous  medium  to  the  gas  (i.e.,  a mixture  of  air  and 
steam)  and  vice  versa.  For  example,  if  the  pores  are  macroscopically  one  quarter 
occupied  with  liquid  water  (degree  of  water  saturation  Sw=  0.25),  then  the  material 
will  effectively  be  less  permeable  to  steam  flow  than  if  it  were  in  a completely 
desaturated  state  (degree  of  liquid  saturation  Sw=  0.0).  Therefore  we  can  no  longer 
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consider  Kg  and  K(  to  be  solely  properties  of  the  solid  matrix,  but  also  we  need  to 

consider  their  dependence  upon  the  extent  to  which  each  fluid  phase  impedes  the  flow 
of  the  other. 

A postulate  commonly  used  to  allow  for  this  interaction  between  the  two  phases 
is  that  the  permeabilities  depend  on  the  fluid  saturations,  the  effective  permeability  of 
each  fluid  phase  increases  with  its  saturation.  At  a macroscopic  scale  we  assume  the 
resulting  permeability  as  a product  of  the  intrinsic  permeability  times  the  relative 
permeability  krx 


W”’  (3.88) 

kf=Kt-km(n,  S ") 

where  kK  stands  for  the  effective  permeability  of  the  solid  matrix  to  a constituent  n at 

the  degree  of  fluid  saturation  S* , n represents  the  porosity  of  the  porous  medium, 
and  the  function  krR  (n,  S*)  are  called  a relative  permeability  function,  which  obeys 

the  bounds  0 < krn  ( n , S" ) < 1 . However,  since  Sg  = 1 - Sw , we  can  define  the  relative 
permeability  as  a function  of  only  the  degree  of  liquid  water  saturation  Sw . Thus, 


K*=KM  s ") 


(3.89) 


Applying  the  new  term  to  both  the  liquid  and  gas  flow  equations,  we  rewrite 
Eqn.  (3.83)  and  Eqn.  (3.85),  for  the  liquid,  as 


rj  v = — 


kK, 


(Vp„-/9wgVz) 


(3.90) 


and  for  the  gas  phase,  as 
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7*v*=_V*  i+±. 

l 


(Vps  -pggVz) 


(3.91) 


For  each  particular  porous  medium,  the  relations  krK(n,  S"  ) are  either  predicted  by 
models  based  on  sorption  equilibrium  or  are  experimentally  determined  in  the 
laboratory.  A curve  of  relative  permeabilities  of  concrete  to  water  and  air  is 
developed  for  this  study  and  will  be  presented  in  Chapter  4. 


3.7.5  Mass  density  of  water 

If  we  set  if/*  = Vs,  i*  = 0 , and  g*  = 0 in  Eqn.  (3.38),  the  mass  conservation  of 
the  liquid  water  phase  can  be  expressed 


Dw(pwvw)  _Q  (3.92) 

Dt 

Assuming  that  pw  is  a function  of  the  fluid  pressure  and  temperature, 
pw  = pw(pw,  T ) , we  can  carry  out  the  differentiation  of  a product  as 


1 Dw(pw)  1 Dw(vw)_  1 (dpw  Dwpw  dpw  DWT' 

~p"  Dt  vw  Dt  ~ pw{  dPw  Dt  + dT  Dt  J 

By  defining  the  bulk  modulus  of  water  Kw  as  the  inverse  of  compressibility 
coefficient  Cw 


(3.93) 
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w 


1 dpw 

PW  dpw 


(3.94) 


and  the  thermal  expansion  coefficient  of  water  aw  as 
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a 


w 


1 dp" 
p"  dT 


(3-95) 


it  follows  that 


pw  Dt  Kw  Dt  w Dt 


1 Dw(pw)  1 Dwpw  . DT 

__  _ 


(3.96) 


Employing  the  equations  (3.94)  and  (3.95)  and  adopting  similar  approach  used  here, 
the  continuity  equation  of  the  fluid  phases  will  be  derived  in  the  later  section. 

3.7.6  The  ideal  gas  law 

Since  the  gas  phase  in  this  study  is  assumed  to  be  a perfect  mixture  of  two  ideal 
gases,  i.e.  dry  air  (superscript  ga)  and  vapor  (superscript  gw),  we  use  the  ideal  gas  law 
that  relates  the  partial  pressure  pa  of  component  a of  the  gas  phase,  the  mass 

concentration  nSgpa  in  Eqn.  (3.57)  and  the  absolute  temperature  0 


where  R represents  the  universal  gas  constant,  Ma  is  the  molar  mass  of  component 
a and  a = ga,  gw . 


(3.97) 


a 


Therefore, 


OR 


(3.98) 


59 


where  the  molar  mass  of  the  moist  air  M is  expressed  with  each  component’s  molar 
mass  and  mass  density 


M 


g 


( cT  1 pga  \ ) 

— — • + — — • 

pg  M pg  M 

\ r lrlgw  H lvlgaJ 


(3.99) 


3.7.7  Fick’s  law  for  diffusion 

At  this  point,  further  physical  interpretation  is  given  for  the  diffusive  mass  flux 
j"  introduced  in  the  macroscopic  mass  balance  law  for  the  gas  phase  of  Eqn.  (3.57). 

Owing  to  the  complicated  system  of  interconnected  passages  (or  channels) 
comprising  the  microstructure  of  porous  media,  variations  in  local  velocity  along 
tortuous  flow  paths  result  from  velocity  distributions  within  each  pore  (both  in 
magnitude  and  direction). 

This  quantity  incorporates  several  types  of  microscopic  effects,  each  of  which 
contributes  to  the  macroscopic  flow  of  a component  (or  a species)  of  a fluid  phase  but 
none  of  which  is  very  amenable  to  laboratory  measurement.  Examples  include 
molecular  diffusion  of  spreading  a solute  within  a fluid  phase,  Taylor  dispersion  of 
longitudinal  spreading  of  the  solute  within  narrow  pore  channels,  stream  splitting  (or 
transverse  spreading)  of  the  solute,  and  the  tortuosity  effects  in  local  velocity  of  the 
fluid  particle  along  tortuous  flow  paths  (Allen  et  al.  1988). 

At  a molecular  scale,  there  is  a substantial  transport  mechanism  in  molecular 
species  dissolved  in  gas  known  as  molecular  diffusion  of  gas,  i.e.  an  additional  mass 
transport  phenomenon,  which  occurs  simultaneously  with  convective  mass  transport 
phenomena  at  the  microscopic  scale  (i.e.  mechanical  dispersion),  is  caused  by 
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molecular  diffusion  resulting  from  variations  in  concentration  of  a gas  species  within 
the  gas  phase  (Bear  and  Bachmat  1998).  Thus,  the  diffusive  flux  includes  two  basic 
transport  phenomena  (1)  molecular  diffusion  and  (2)  mechanical  dispersion. 


Je  Jmole  Jmech 


(3.100) 


where  jmole  and  jmech  represent  the  molecular  diffusion  and  mechanical  dispersion 
fluxes,  respectively. 

Actually,  the  separation  between  the  two  processes  is  rather  artificial,  as  the 
transport  mechanism  of  the  gaseous  phase  includes  both  processes  in  an  inseparable 
form  named  hydrodynamic  dispersion.  Thus,  the  diffusive  mass  flux  j“  in  Eqn. 
(3.57)  can  be  written  in  a form  combining  the  two  processes 


j“  = -nSgpgr Dg  grad 

where  Dg  is  the  effective  dispersion-diffusion  coefficient  tensor,  r is  the  tortuosity 

factor,  and  a is  the  diffusing  component  in  the  gaseous  phase  . This  is  known  as 
Fick ’s  law  of  binary  diffusion.  However,  since 


(3.101) 
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then. 
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(3.103) 


Or,  by  using  the  ideal  gas  law,  it  can  be  also  expressed 
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(3.104) 


where  a represents  component  of  the  gaseous  phase,  a - ga,  gw,  which  stands  for 
dry  air  and  vapor,  respectively.  It  should  be  noticed  that  the  gas  diffusion  is  an 
additional  mass  flux  due  to  a mass  concentration  contributable  to  the  gas  phase  flow 
even  in  the  absence  of  a pressure  gradient. 


3.7.8  Fourier’s  law 

A constitutive  equation  used  for  the  heat  flux  by  conduction  is  the  generalized 
Fourier’s  law  that  is 


q = §rad  9 


(3.105) 


where  q is  the  averaged  heat  flux  of  the  multiphase  medium  that  is  the  sum  of  the 
partial  heat  fluxes  qK  in  a phase  n , and  Keff  is  the  effective  thermal  conductivity 


tensor  for  the  multiphase  medium 


*tf=2X**  K = 5>  wand  g 

n 


(3.106) 


where  s,  w and  g represent  the  solid,  the  liquid  water,  and  the  gas  phase, 
respectively. 


3.7.9  Sorption  isotherms 

Microstructure  characteristics  of  a porous  medium  have  been  expressed  at  the 
macroscopic  level  by  transfer  coefficients  such  as  permeability,  diffusivity  and 
conductivity,  which  have  direct  influence  on  mechanical  behavior.  In  particular, 
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textural  parameters  such  as  the  porosity  and  the  pore  size  distribution  are  often 
thought  to  be  crucial  in  the  movement  of  liquid  and  gaseous  water  in  the  pore 
network  when  the  porous  medium  is  exposed  to  natural  drying  conditions. 

Determination  of  the  pore  size  distribution  and  thus  the  porosity  of  the  porous 
medium  can  be  achieved  from  a laboratory  test  known  as  water  vapor  desorption  and 
adsorption  test  (e.g.  if  an  oven-dry  porous  sample  is  exposed  to  moist  air,  the  weight 
of  the  sample  increases  due  to  the  moisture  adsorption  on  the  inner  surfaces  of  the 
pores).  Sorption  isotherms  can  be  used  to  represent,  (1)  the  actual  water  content  under 
natural  drying  conditions  which  depend  on  the  initial  boundary  conditions  in  the 
numerical  modeling,  and  (2)  the  relative  permeability  that  is  a function  of  the  porosity 
and  the  degree  of  water  saturation. 

The  postulate  of  sorption  isotherms  can  be  summarized  as  follows  (1)  the 
capillary  pressure  pc  is  defined  as  the  pressure  difference  between  the  gaseous  phase 
(a  mixture  of  vapor  and  dry  air)  and  the  liquid  phase  (water)  as 


PC=Pg-P 


(3.107) 


(2)  a relationship  between  the  relative  humidity  (RH)  and  the  degree  of  water 
saturation  Sw  can  be  determined  by  water  sorption  isotherms  (i.e.  an  equilibrium 
curve  of  water  content  versus  RH)  in  a laboratory  test,  and  (3)  the  degree  of  water 
saturation  Sw  are  assumed  to  be  a function  of  the  capillary  pressure  and  the 

temperature  difference  T above  a reference  temperature  0{) 


Sw  = Sw(pc  ,T) 


(3.108) 
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where  T -0-0o. 

3.7.10  Total  stress  and  effective  stress  in  the  solid  phase 

Combining  the  total  stress  tensor  defined  in  Eqn.  (3.71)  and  the  stress  tensor 
written  in  Eqn.  (3.79),  we  rewrite  the  total  stress  tensor  as 


a'°'al  = Tjscrs  + + rjh 7*  = (1  - n)as  - n{SwPw  1+  SgPg  I ) 


(3.109) 


or 


<J  = TJS(7S  + T]waw  + Tjg(Tg  - (1  - n)as  - n(Ps I ) 


(3.110) 


where  <7lolal  is  the  total  stress  tensor  and  Ps\  is  the  mean  hydrostatic  pressure  tensor 


of  both  the  fluids  surrounding  the  grains,  that  is 


Ps=SwPw  + SgPg 


(3.111) 


By  rearranging  the  terms,  the  constitutive  law  of  the  solid  phase  can  be  expressed 
through  the  concept  of  Terzaghi’s  effective  stress  creff 


cr,olal  = (1  - n)as  + (1  - n\Ps  I ) - n(ps  I ) - (1  - n){Psl ) 
= (\-n)(as  +Psl)-Psl 

= (7eff  - Psl 


(3.112) 


Therefore  the  effective  stress  that  is  responsible  for  the  deformation  in  the  solid  phase 


is 


a‘s  = a"M“  + pf 


(3.113) 


where  aeS  is  known  as  Terzaghi’s  effective  stress. 
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3.8  Physical  Approximations  and  Field  Equations 

Having  introduced  the  constitutive  equations,  we  now  transform  the 
macroscopic  balance  laws  to  obtain  the  field  equations,  which  will  be  solved 
numerically  in  the  following  chapters.  Although  each  term  in  the  general  macroscopic 
equations  constitutes  a mathematical  model  that  describes  the  contribution  of  a 
certain  phenomenon  to  the  process  considered  herein  (i.e.,  thermo-mechanical 
behavior  of  concrete  under  multi-phase  coupled  heat  and  mass  transport),  only  a 
small  number  of  phenomena  may  turn  out  to  be  dominant  in  comparison  to  others. 

In  what  follows,  we  will  discuss  the  key  assumptions  that  need  to  be  made  for 
deletion  of  terms  representing  non-dominant  effects  in  the  considered  process.  This 
allows  us  to  produce  simplified  balance  equations  that  may  be  rather  less  complicated 
yet  more  manageable  to  solve  in  the  case  of  non-isothermal  two-phase  flow  in 
deforming  porous  media,  with  relatively  little  effect  on  its  solution,  i.e.,  within  a 
permissible  range  of  error. 

3.8.1  Continuity  equation  for  the  fluid  phases 

By  assuming  that  the  solid  velocity  is  negligible  in  the  macroscopic  mass 
balance  equations  for  the  solid  phase,  Eqn.  (3.45)  becomes 

(1  -n)dps  dn  0 (3.114) 

ps  dt  dt 

For  the  liquid  water,  Eqn.  (3.48)  can  be  rewritten  in  terms  of  the  material 
derivative  with  respect  to  the  solid  by  using  the  definition  of  material  derivatives 
(refer  to  Eqn.  (3.36))  and  the  vector  identity  (also  see  Eqn.  (3.43)) 
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D\nSwpw) 

Dt 


+ v™ -grad  ( nSwpw ) + nSwpw div (v5  + vws)  = -m 


(3.115) 


However,  since  v'  = 0 , the  relative  velocity  of  the  liquid  water  to  the  solid  is 
replaced  by  the  phase’  velocity.  This  leads  to  rewriting  the  mass  balance  equation  of 
the  liquid  water  in  the  form  of 


D\nSwp”) 

Dt 


+ vw«grad  (nS  w pw  ) + nSwpw  div  vw  =~m 


By  applying  the  vector  identity,  it  becomes 


(3.116) 


D\nSwp") 

Dt 


+ div  ( nSwpwvw ) = -m 


(3.117) 


Small  displacements  are  assumed  for  the  solid  phase.  Thus,  the  material  time 
derivative  Ds(»)/  Dt , is  replaced  by  a partial  time  derivative  d(»)/dt 


d{~nSwP  ) + div ( nSwpwvw ) = -m 
dt 


(3.118) 


Expanding  and  dividing  both  sides  of  the  equation  by  Swpw , we  find 


dn  n dS  n dpw  1 „ w w.  m 

dt  Sw  dt  p"  dt  Swpw  Swpw 


(3.119) 


Summing  Eqn.  (3.1 19)  to  Eqn.  (3.114),  we  have 


Ps  dt  sw  dt  p"  dt  swP * 


PWS» 


(3.120) 


This  is  the  continuity  equation  for  liquid  water. 
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By  referring  to  Eqn.  (3.94)  and  Eqn.  (3.95)  and  defining  the  thermal  expansion 
coefficient  of  the  solid  skeleton  a . as 


= -■ 


1 dp5 
ps  dT 


(3.121) 


the  continuity  equation  of  liquid  water  can  be  rewritten 


(1  ~n) 


( dT\ 
-a, — 
y dtj 


n dSw 

+ - + n 

SL  dt 


r 1 dpw  dT ' 

— — 

vkw  dt  dt 


(3.122) 


+ — 5— di \(nSwpw  vw)  = — 

SWP  P sw 


Rearranging  the  terms  and  then  multiplying  Sw  on  the  both  sides,  we  have 


Sw(-(l-n)as -nccw)—+  nS"  dPw  vw)  - 

V 5 ’ dt  Kw  dt  dt  pw  w 


m 

v 


(3.123) 


or 


dT  nS  dp  dS  1 ,.  w w m 

or — + — -~^  + n — - + — di \(nSwp  v ) = 

dt  Kw  dt  dt  pw  pw 


(3.124) 


where  asw  is  defined  as  the  relative  thermal  expansion  coefficient  of  the  solid  to  the 
liquid  water 


=Sw(-(\-n)as-naw) 


(3.125) 


Finally,  by  substituting  the  generalized  Darcy’s  law  from  Eqn.  (3.90)  into  the 
continuity  equation  of  the  liquid  water,  it  follows  that 


67 


dT  _ nSw  dpw  t dSw  1 
a„,  — + — -~L^-  + n — - + — div 


dt  K dt 


dt 


p'^C-gradp.  + p'gVz) 


(3.126) 

Pw 


Similarly,  the  continuity  equation  of  the  gaseous  phase,  i.e.  moist  air,  can  be 
derived  in  a same  manner  as  that  used  to  derive  the  mass  balance  of  the  liquid  water, 
with  the  result  being 


(1  ~n) 
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dt  . 
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S'  dt 
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P'S, 


However,  since 


^ = _aSs,  (3.128) 

dt  dt 


we  have 
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g dt 
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-di v(nS  p8  Vs)  — 


rh  (3.129) 
P8Sg 


This  equation  can  also  be  split  into  two  parts  due  to  the  species  mass  balance  law 
defined  in  Eqns.  (3.55)  and  (3.57).  For  dry  air,  Eqn.  (3.129)  is  rewritten 


(1-h) 
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Q> 
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Sg  dt 
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(3.130) 


For  vapor,  Eqn.  (3.129)  is  rewritten 
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(l-») 
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8T_ 

dt 


n dSw  n d f PgvMl ^ 


SK  dt  p8”  dt\  OR 


+ —1  ^ [divjf  +div(«^P^  vg)]  = 


(3.131) 


m 


S*P 


Substituting  Fick’s  law,  Eqn.  (3.104),  and  the  generalized  Darcy’s  law,  Eqn. 
(3.85),  into  the  last  two  continuity  equations  for  the  components  of  the  gaseous  phase, 
we  have  the  following  expressions 


(l-») 
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n dSw  n d 
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(3.132) 
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and  for  vapor, 
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(3.133) 
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where,  again,  the  first  three  terms  on  the  left  side  account  for  the  accumulation  of  the 
components  known  as  the  mass  concentration  in  the  gaseous  phase,  i.e.  moist  air,  the 
forth  term  models  the  diffusion  of  each  component  with  respect  to  the  mean 
macroscopic  flow  of  moist  air  in  the  porous  medium,  the  last  term  represents 
advection  of  the  components  with  the  mean  macroscopic  flow  velocity,  and  the  right 


69 

side  accounts  for  the  exchange  of  vapor  to  the  liquid  water  phase  in  the  porous 
medium. 

3.8.2  Enthalpy  balance  equation  for  the  multiphase  medium 

Recall  the  macroscopic  energy  balance  presented  in  Eqn.  (3.76), 


f pn  ^ — = rfp*R*  - 77*  V-q'  + r/ ■ Vv*  + rf  pK  (h*  + g. v*  ) 


(3.134) 


Although  the  macroscopic  energy  balance  law  provides  a good  basis  for 
understanding  the  nature  of  heat  transfer  in  a porous  medium,  we  need  to  develop  a 
more  manageable  field  equation  to  solve  numerically. 

For  further  development,  the  internal  energy  of  a phase  may  be  used  as  a state 
variable,  that  is,  a variable  that  depends  only  on  the  state  of  the  phase  and  not  on  any 
process  that  produced  that  state.  Therefore,  we  are  free  to  define  additional  state 
variables  that  are  combinations  of  existing  state  variables.  These  new  variables  will 
often  make  the  analysis  of  a system  much  simpler.  A useful  additional  state  variable 
is  the  specific  enthalpy  (H)  that  is  defined  to  be  the  sum  of  the  specific  internal  energy 
(. E)  plus  the  product  of  the  pressure  ip)  and  the  specific  volume  ( 1 Ip).  Then  we 
express  the  energy  balance  by  means  of  the  specific  enthalpy  H defined  as 


(3.135) 


where  E * is  the  intrinsic  specific  energy  of  a control  volume  of  a phase  n and  — is 


the  intrinsic  specific  volume.  The  differential  of  the  specific  enthalpy  is  then 
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dH*  = dE*  + pKd 


kP"  j 


1 , 

+ -JTdp« 
P 


(3.136) 


By  rearranging  the  terms,  it  follows  that 


f i \ 


dEn  = dET  -pKd 


P ) P 


—JdP* 


(3.137) 


Therefore,  the  left  hand  side  of  Eqn.  (3.134)  can  be  written 


„ XD*E'  , XD*H*  , 1 D*p*  „ Dnp 

v p —^—=n  P —p—+Ti  P*——£r — n 
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Dt 
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pn  Dt 

— V 
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Dt 


(3.138) 


Let  us  define  the  terms  ( a ) and  ( b ) of  the  right  hand  side  of  the  equation  as  the 
following.  Since  the  specific  enthalpy  is  a function  of  the  absolute  temperature  and 
the  pressure 


H*  =H*(0,  pK) 


(3.139) 


its  differential  can  be  written  as 


dH * = 


v 36*  j 


d9*  + 


( dH*^ 


V dp*  Jo 


dP, 


(3.140) 


By  using  the  definition  of  the  specific  heat  at  constant  pressure 


c;  = 


f dH*^ 


(3.141) 


and  the  following  thermodynamic  identity  (Hudson  1 996) 
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(3.142) 


Therefore  the  term  (a)  of  Eqn.  (3.138)  becomes 


---  = rK 
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(3.143) 


p. 


Dt 


For  the  material  time  derivative  of  the  phase  intrinsic  density,  we  expand  and 
divide  Eqn.  (3.48)  by  rf  p* 


1 Vp*  . , 1 DKrf 

■ = e -div  v — 


(3.144) 


pK  Dt  rf  Dt 

Substituting  Eqns.  (3.143)  and  (3.144)  into  Eqn.  (3.138),  we  have  the  left  hand 
side  of  Eqn.  (3.134)  as  the  following 


r\7T  V)n  f)K 

rf  p*  — = rf  pnCnp  — - — rfQn 
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(3.145) 


-TfP*e*-rfP, div 
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Next,  we  consider  the  mechanical  energy  term  aK  : Vv*  of  Eqn.  (3.134)  as  the 
sum  of  the  deviatoric  part  and  the  hydrostatic  part 


cr*  : Vv*  = Tn  : Vv*  - /^V-v* 


(3.146) 


If  we  set  the  exchange  of  internal  energy  (e.g.  latent  heat)  between  phases  due 
to  phase  change  p*  R*  of  Eqn.  (3.134)  as 
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-71  T)K  _/T  _ n T?n 

p K = p e t 


(3.147) 


then  the  latent  heat  pneK  HK  of  evaporation  for  the  fluid  phases  are  defined  by  using 
the  definition  of  the  specific  enthalpy  of  Eqn.  (3.135)  as 


n n j j k n „7i  j~>tt  , n 

pen  = p e t + p^e 


(3.148) 


By  substituting  the  latent  heat  term  and  the  mechanical  energy  term  into  Eqn.  (3.134), 
the  macroscopic  enthalpy  equation  is  formed 
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where  (a)  represents  the  rate  of  change  of  heat  stored  in  the  phase,  ( b ) is  the  heat 
source,  (c)  represents  the  internal  heat  flux  in  the  phase,  (d)  accounts  for  the  latent 
heat  of  evaporation  or  condensation  for  the  fluid  phases,  (e)  mechanical  work  by 
volume  fraction  change,  (/)  represents  external  work  done  by  body  forces,  (g)  is 
viscous  dissipation  by  the  deviatoric  stresses,  and  ( h ) accounts  for  mechanical  work 
due  to  density  variation  caused  by  temperature  changes. 

For  further  development,  we  assume  that  a state  of  local  thermodynamic 
equilibrium  in  all  the  phases  exists.  This  means  that  the  temperatures  of  the  phases 
are  equal  at  each  point  (or  REV)  in  the  multi-phase  system 
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es  =ow  = e8 


(3.150) 


Also,  the  terms  related  to  viscous  dissipation  and  mechanical  work  done  by  the 
density  variation  and  the  volume  fraction  change  are  assumed  to  be  negligible. 

Recalling  Fourier’s  law  and  the  definition  of  material  time  derivative,  we  find 
the  macroscopic  enthalpy  equation 


dT 


(p  C,)tf  ~-+(nS,p’C;w’  + nSgp‘C‘v‘ygad  T 


(3.151) 


div  (Keff  grad  T)  = -mAFI, 


vapor 


where 


(p  Cp)eff  = (1  -n)psCsp+nSwp"C;  + nSgp*Cl 

Keff  = 0 " «)<  + nS^wp  + nSgKl  (3.152) 

A Hmpor=H*"-H” 

The  temperature  difference,  T = 0 - 0O , is  substituted  for  the  absolute  temperature  0 
in  all  derivatives  where  the  reference  temperature  0O  is  fixed. 

A physical  interpretation  of  Eqn.  (3.151)  can  be  given  that  the  first  term  on  the 
left  represents  rate  change  of  heat  stored  in  the  multiphase  medium,  the  second  term 
accounts  for  heat  flow  by  convection,  the  third  term  represents  the  internal  heat  flux 
by  conduction,  and  the  term  on  the  right  hand  side  represents  the  enthalpy  change  due 
to  evaporation  of  the  liquid  water. 

3.9  Initial  and  Boundary  Conditions 

The  initial  and  boundary  conditions  for  a partially  saturated  porous  medium 
subjected  to  high  temperature  are  now  considered.  The  initial  conditions  describe  the 
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field  variables,  i.e.  displacements,  water  pressures,  gas  pressures,  degrees  of  water 
saturation  (or  degrees  of  gas  saturation),  and  temperatures  at  time  zero 

w = uo  Pw  = pI  Pg=P°g  Sw=sl  T = T0  in  Q and  on  <9Q  (3.153) 

where  heat  and  mass  transfer  coupled  with  mechanical  displacements  occurs  in  a 
partially  saturated  porous  material  occupying  the  volume  Q with  a closed  surface 
boundary  dQ  over  a time  interval  / = [0, 7’] 

We  introduce  the  Dirichlet  type  of  boundary  conditions,  i.e. 


u-u  on  <3Q], 

Pw=K  on 
pg=pg  on  5Q' 

T = f on  8D!t 


(3.154) 


For  the  Neumann  type,  first  the  traction  t boundary  condition  for  the  mechanical 
stresses  am  can  be  written 


<7m  • n = t on  dQ2u 


(3.155) 


where  n is  a unit  normal  vector  to  the  boundary  surface  and  the  respective  flux 
boundary  conditions  of  each  phase,  i.e.  for  the  water  flux 
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• n = qw  + qgw  ondQ2w  (3‘156) 


ypg) 


75 


for  the  gas  flux 


• n = qga  on  SQ2  (3‘157) 


and  for  the  heat  flux 
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where  dQ  = 3Q1  + 3Q2  and  q is  the  imposed  mass  (or  heat ) flux  normal  to  the 
boundary.  These  equations  are  the  natural  boundary  conditions  for  the  macroscopic 
balance  equations,  the  solution  to  which  is  obtained  via  a ‘weak’  formulation  (e.g. 
using  the  weighted  residual  methods)  of  the  problem  as  is  done  in  Chapter  6. 

In  heat  flux  boundary  conditions,  three  basic  mechanisms  are  involved  in  heat 
transfer  between  boundaries  from  one  body  to  another  conductive,  convective,  and 
radiant  heating.  If  we  attach  a metal  plate  to  another  hot  plate,  we  eventually  notice 
that  the  attached  plate  becomes  warmer.  This  form  of  heat  transfer  is  called 
conduction.  In  a warm  fluid  flowing  in  a pipe,  heat  of  the  fluid  is  transferred  by  the 
mechanism  called  convection  to  the  pipe.  In  contrast  to  the  mechanisms  of  conduction 
and  convection,  where  heat  transfer  through  a medium  is  involved,  heat  may  also  be 
transferred  by  electromagnetic  waves  that  are  propagated  as  a result  of  a temperature 
difference.  This  is  called  radiation. 

Radiation  is  a significant  mode  of  heat  transfer  in  building  fire  environments,  it 
is  particularly  important  to  modeling  enclosure  fires.  Besides,  radiation  tends  to 
dominate  the  other  two  modes  of  heat  transfer  in  room  environments  where  high 
temperatures  are  attained  with  fires  or  hot  smoke  layers.  Therefore,  modeling  of 


76 

radiation  boundary  conditions  between  a heat  source  and  a body  is  of  particular 
interest. 

In  a building  fire  scenario,  the  conductive  heat  flux  boundary  condition  is 
considered  in  the  mathematical  expression  of  heat  flux  boundary  condition  in  Eqn. 
(3.158).  Two  basic  numerical  approaches  can  be  considered  for  modeling  the  radiant 
boundary  conditions  (Consolazio  et  al.  1998)  (1)  computing  heat  fluxes  in  boundary 
elements  directly  from  the  heat  source,  and  (2)  converting  the  radiant  heat  fluxes  on 
the  boundaries  into  a conductive  boundary  condition.  The  latter  approach  is  adopted 
in  this  research.  Thus,  in  the  mathematical  expression,  only  conductive  boundary 
condition  is  included  hereby.  Additional  details  describing  the  second  approach  for 
computing  radiative  heat  transfer  between  the  bounding  surfaces  of  concrete  and 
point  source  fires  are  obtained  from  Patankar  (1980). 


CHAPTER  4 

NUMERICAL  MODELING  OF  TRANSPORT  PHENOMENA  IN  REINFORCED 
CONCRETE  EXPOSED  TO  ELEVATED  TEMPERATURES 

The  objective  of  this  chapter  is  to  develop  a finite  difference  model  that  simulates 
coupled  heat  and  mass  transport  phenomena  in  reinforced  concrete  exposed  to  rapid 
heating  conditions  such  as  fires.  A mathematical  and  computational  model  to  simulate  the 
multi-dimensional  thermo-hydrological  response  of  reinforced  concrete  structural 
elements  is  presented.  Key  parameters  that  describe  thermo-hydrological  material 
behavior  of  concrete  are  also  presented.  Effects  of  steel  reinforcements  on  coupled  heat 
and  mass  transfer,  which  can  be  potentially  associated  with  concrete  spalling,  are 
discussed.  The  domain  and  time  history  distributions  of  temperature,  pore  pressure,  and 
degree  of  saturation  are  illustrated  as  predicted  under  various  thermal-loading  conditions. 

4.1  Overview 

The  lower  permeability  of  HSC  limits  moisture  movement  in  the  heated  concrete 
and  results  in  the  development  of  large  internal  pore  pressures  (Kodur  1997,  Consolazio 
et  al.  1998,  Kodur  1999).  In  addition,  material  properties  (e.g.,  material  strength  and 
thermal  expansivity)  and  mechanical  behavior  (e.g.,  stress  and  strain  fields)  can  be 
significantly  influenced  by  the  thermodynamic  state — as  characterized  by  variables  such 
as  pore-pressure,  temperature,  and  saturation — of  reinforced  concrete  elements. 

Thus,  the  goal  in  this  chapter  is  to  quantify  thermodynamic  state  variables  for 
concrete  under  a variety  of  different  heating  conditions.  By  quantifying  state 
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variables,  greater  insight  into  the  causes  of  phenomena  such  as  explosive  spalling  can 
be  obtained.  Focus  in  this  chapter  is  placed  primarily  on  transient  heat  and  mass 
transfer  occurring  inside  heated  concrete  rather  than  on  prediction  of  resulting 
stresses — the  topic  of  later  chapters.  In  addition  to  insight  into  flow  characteristics  of 
HSC,  the  effects  of  steel  reinforcement  on  hydro-thermal  behavior  of  concrete  are 
also  investigated  in  order  to  quantitatively  identify  the  influence  of  heterogeneity  on 
the  thermodynamic  state  of  the  system.  For  example,  steel  reinforcement  can  create 
impermeable  regions  that  alter  the  migration  of  moisture  and  water  vapor  through  the 
heated  concrete  matrix.  An  analytical  model  including  key  flow  parameters  for  high- 
strength  concrete  that  simulates  transient  fluids  (i.e.,  water  vapor  and  liquid  water) 
and  heat  flow  is  therefore  developed  in  this  chapter. 

Formulating  realistic  mathematical  representations  of  dominant  physical 
phenomena  is  a primary  component  of  the  model  development  process  described 
here.  The  TOUGH2  (Pruess  1991)  code — a general  purpose  finite  difference 
numerical  simulator  for  multi-phase  fluid  and  heat  flow — is  used  as  the  basis  for 
numerical  modeling.  The  mathematical  formulation  implemented  in  TOUGH2 
consists  of  mass  and  thermal  energy  balance  equations  that  are  used  to  solve  for  three 
primary  variables  pore  pressure,  degree  of  gas  saturation,  and  temperature.  However, 
in  this  study  modifications  are  made  to  the  code  so  that  important  flow  characteristics 
specific  to  concrete  behavior  are  properly  taken  into  consideration.  For  example, 
relative  permeability  functions  and  slip  flow  relationships  are  developed  specifically 
for  cement-based  materials  and  are  implemented  into  TOUGH2.  Implementation  of 
these  new  features  is  necessary  in  order  to  properly  describe  the  flow  of  the  various 
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fluid  phases  that  migrate  inside  heated  concrete.  The  proposed  model  is  capable  of 
accounting  for  the  simultaneous  movement  of  both  the  gaseous  and  liquid  phases, 
their  transport  of  heat,  and  phase  transitions  between  liquid  water  and  gaseous  water 
vapor. 


Having  introduced  the  mass-  and  energy  balance  equations  in  generic  form  in 
Chapter  3,  we  now  review  the  governing  equations  of  coupled  heat  and  mass 
transport  to  be  solved  for  three  primary  variables  pressure,  temperature,  and 
saturation.  The  mass-  and  energy-balance  equations  can  be  expressed  in  the  integral 
form  of 


where  part-a  represents  the  rate  of  change  of  mass  accumulation  for  fluid  phases  (if 
y- 1 or  2)  or  the  rate  of  change  of  heat  accumulation  for  the  multiphase  continuum 
( y = 3 ),  part-6  represents  the  mass  flux  of  a fluid  phase  ( y =1  or  2)  or  the  heat  flux 
( y = 3 ) into  U through  surface  area  dU  , and  part-c  represents  the  external  supply  of 
energy  attributable  to  gravitational  effects. 

More  specifically,  the  mass  accumulation  term  M(r)  in  Eqn.  (4.1)  can  be 
expressed 


4.2  Governing  Equations 


(4.1) 


(y  = 1 : water;  y -2:  air;  y = 3 : heat) 
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where  the  water  specie  can  be  written  in  terms  of  liquid  and  gas  phases  as 


That  is,  the  mass  accumulation  of  water  species  is  equal  to  the  sum  of  masses  of 
liquid  water  and  water  vapor.  For  the  air  specie, 


where  n is  the  total  porosity,  S„  is  the  degree  of  saturation  in  each  phase,  pn  is  the 
intrinsic  phase  density,  and  represents  the  mass  fraction  of  component  y in 
phase  n , which  can  be  written 

— P (4 

* p* 

The  heat  accumulation  term  accounts  for  thermal  energy  stored  in  the 
multiphase  continuum. 


where  C'p  is  the  specific  heat  of  the  porous  matrix  and  En  represents  the  specific 
internal  energy  of  each  water  phase  n ( n - w,  g ) 

The  mass  flux  term,  F(r) , in  Eqn.  (4. 1)  is  a sum  of  the  mass  fluxes  summed  over 


(4.3) 


(4.4) 


M«>=(1  -n)p-CrT  + n £ S,p’E, 


(4.6) 


each  phase 
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F<r)  _ ^ F00  (4.7) 

n=w,g 

Based  on  a generalized  Darcy’s  law  modified  for  multi-fluid  flow  and  slip  flow, 
the  mass  flux  of  water  species  ( y = 1 ) can  be  expressed  as 


Fd)  _ Fd)  , Fd) 

x w A g 

= %P"X™(VA-^' gVz)-K 

M Ms 


f h ^ 

1 + -^ 

V P'J 


(4.8) 


(Vp~pggVz) 


and  the  mass  flux  of  dry  air  ( y = 2 ) can  be  written  in  the  form 


F(2)  = -K^-pgX?) 


p 

g y g 


f h ^ 

1 + A 

PtJ 


(Vpg-  psgVz)~  p8zT)g  grad  (x(g2)) 


(4.9) 


where  K is  the  intrinsic  permeability,  krtt  is  the  relative  permeability  to  phase  n , pK 

is  the  dynamic  viscosity,  z is  the  tortuosity  factor,  and  Dg  is  the  effective 

dispersion-diffusion  coefficient  for  diffusing  dry  air  in  the  gaseous  water  phase  (i.e., 
Pick ’s  law  of  binary  diffusion). 

The  total  heat  flux  term,  F(3),  in  Eqn.  (4.1)  also  consists  of  two  parts,  heat 
conduction  through  the  multiphase  continuum  as  a whole  and  heat  convection  by 
mass  fluxes  of  the  fluid  phases  (liquid  water  and  moist  air) 


F,!>=-v^  grad  T+£ //">!*/>  (4.10) 

r-l.l 

1,2 


where  Hf]  represents  the  specific  enthalpy  of  component  y in  phase  n , Kmg  is  the 


volume-averaged  thermal  conductivity  tensor  for  the  multiphase  medium 
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Karg  = Z ^ ;r  = s,  w and  g (4.1 1) 

n 

where  rf  is  the  volume  fraction  of  phase  n , and  kk  is  the  thermal  conductivity 
tensor  of  phase  n . 

4.3  Discretization  of  Integrated  Finite  Difference  Method 

The  integral  mass-  and  energy-balance  equations  in  Eqn.  (4.1)  are  solved 
directly  from  the  integral  formulation  known  as  the  ‘ integrated  finite  difference 
method"  (IFDM)  (Narasimhan  and  Witherspoon  1976).  Consider  two  volume 
elements  V(  and  Vm  with  the  interface  surface  segment  A(m  where  the  distances 
between  the  centroids  of  the  volumes  and  the  interface  are  denoted  by  Df  and  Dm  , 
shown  as  Figure  4-1. 


Figure  4-1  A schematic  sketch  of  spatial  discretization 
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The  mass  flux  of  the  water  components  ( y = 1 ) can  be  expressed  in  discretized 
form  as 


F<»  - 


-K 


wjm 


1 

E 

1 

1 

l ^ ) 

fm 

L J 

(4.12) 
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fm 
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s 

1  

where  the  subscripts  ( fm ) indicate  the  averaged  value  at  the  interface  between 
volume  elements  V(  and  Vm  and  is  based  on  the  macroscopically  volume-averaged 

quantities  within  V(  and  Vm , KeK  (m  represents  the  weighted  value  of  the  intrinsic 

permeability  to  phase  n when  the  two  volume  elements  are  composed  of  different 
materials  or  composed  of  the  same  material  but  which  are  partially  saturated  to 
different  degrees  for  phase  n , pKl  represents  the  average  pressure  of  phase  n in 

volume  element  Vf , and  D(m  represents  the  distance  between  the  centroids  of  volume 
elements  V(  and  Vm . Gravitational  effects  are  omitted  because  the  focus  here  is  on 
localized  behavior  (e.g.  spalling)  and  gravitational  effects  are  not  of  significant 
importance  in  such  situations. 

The  effective  permeability  to  water  Kewlm  at  the  interface  between  volume 
elements  V(  and  Vm  is  defined  as  (Narasimhan  and  Witherspoon  1976) 


KWJDm+Kw<mDe 


(4.13) 


whereas  the  effective  permeability  to  moist  air  is  defined  as 
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(4.14) 


Kg,eDm  + KgjnDC 


where  gas  permeability  within  each  volume  element,  including  Klinkenberg’s  effect, 
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Mass  flux  of  the  dry  air  can  be  expressed  in  discretized  form  as 


(4.15) 
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(4.16) 


where  pfm  is  an  averaged  mass  density  of  the  gaseous  phase  in  volume  elements  V( 
and  Vm , r(m  is  an  average  value  of  the  tortuosity  factor,  Ds  (m  is  an  average  value  of 
the  effective  dispersion-diffusion  coefficient,  and  X{*\  represents  the  mass  fraction  of 
component  y in  phase  n within  volume  element  V( . 

Similarly,  for  the  heat  flux  term  F<3)  in  Eqn.  (4.1),  the  heat  flux  can  be 
discretized  as 


F(3)  = -Ke 

fm  avgff  m 


Tt-Tm 

Dftn 


*=1,2 


(4.17) 


where  Kemg  tm  is  the  harmonic  mean  (Patankar  1 980)  of  conductivity  coefficients 
Kmg,i  an(l  Kavg, m > and  is  denoted  the  “effective”  conductivity,  Tt  and  Tm  are 
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temperatures  in  volume  elements  Vf:  and  Vm  at  a local  thermodynamic  equilibrium 
state,  and  Hn  represents  an  averaged  value  of  the  sum  of  specific  enthalpies  of  all 
the  components  (or  species)  y in  phase  n . 

4.4  Key  Parameters  of  Fluid  and  Heat  Flow  in  Concrete 

In  modeling  simultaneous  flow  of  heat  and  fluids  through  concrete,  the  effects 
of  moisture  movement  on  key  flow  parameters  such  as  permeability  and  conductivity 
must  be  included.  For  example,  permeability  to  moist  air  in  the  transport  process  is 
affected  by  presence  of  liquid  water  in  pores.  Also,  thermal  conductivity  of  the 
concrete  is  affected  by  degree  of  water  saturation  during  the  transient  hydro-thermal 
process.  While  TOUGH2  (Preuss  1991)  is  an  excellent  tool  for  the  numerical  analysis 
of  the  continuum  phenomena  of  mass  and  heat  transfer  in  concrete,  reliable  data  for 
key  parameters  used  in  the  constitutive  models  (e.g.  flow  properties  of  a particular 
phase)  are  still  needed  to  make  valid  predictions  by  simulation.  This  study  makes  use 
of  hydro-thermal  properties  collected  from  field  tests  reported  in  literature. 

4.4.1  Relative  permeability  of  concrete 

The  model  presented  here  accounts  for  the  involvement  of  fluid  phases  in  heat 
transfer  that  includes,  (1)  flow  of  both  liquid  and  gas  being  accompanied  by  heat 
convection,  and  (2)  evaporation  of  free  water  content  associated  with  latent  heat. 
Temperature  data  obtained  from  analyses  including  these  effects  are  thought  to  be 
more  accurate  than  similar  data  predicted  by  traditional  thermal  analysis  models 
which  often  do  not  account  for  the  changes  of  heat  conductivity  introduced  by 
changes  in  degrees  of  water  saturation. 
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We  define  the  permeability  as  an  intrinsic  permeability  (also  called  absolute 
permeability)  for  use  in  numerical  simulations.  Although  rigorous  experiments  by 
many  researchers  have  been  conducted  to  measure  the  permeability  of  concrete  to 
either  gas  or  liquid  (Loughborough  1966,  Hope  and  Malhotra  1984,  Nagataki  et  al. 
1985,  McVay  and  Rish  1995,  Reinhardt  and  Dinku  1996,  Sugiyama  et  al.  1996),  only 
a limited  number  of  studies  (Bamforth  1987a  and  b,  Whiting  1988,  Dhir  et  al.  1989) 
were  found  to  provide  data  for  both  liquid  and  gas  permeability  tested  with  the  same 
concrete  mixtures.  Additionally,  consistency  in  preparing  the  samples  (in  particular, 
curing  procedures  and  periods  of  time  are  important  in  air  permeability  measurement 
that  is  substantially  reduced  by  moisture  content  (Reinhardt  and  Dinku  1 996, 
Nagataki  et  al.  1985,  Hope  and  Malhotra  1984))  with  a similar  composition  of  the 
mixture  including  W/(C+P)  and  Agg/(C+P)  was  the  primary  consideration  in 
collecting  data  for  the  key  parameters.  Therefore,  the  parameters  used  in  this  study 
may  not  characterize  all  possible  concrete  types,  but  they  can  be  considered  to  be 
representative  of  normal-weight  HSC  often  found  in  practice. 

Because  the  macroscopic  permeability  of  concrete  is  a function  of  both  the  solid 
matrix  and  the  particular  fluid  used  to  experimentally  measure  the  permeability 
(Klinkenberg  1941,  Nagataki  et  al.  1985,  Dhir  et  al.  1989,  McVay  and  Rish  1995), 
neither  gas  permeability  nor  liquid  permeability  is  believed  to  be  the  single  absolute 
measure  of  the  intrinsic  permeability  of  the  solid  skeleton  of  concrete.  Thus,  for  a 
complete  description,  fluid  flows  in  concrete  will  be  characterized  by  three 
parameters  the  intrinsic  gas  permeability  Kg  , the  intrinsic  water  permeability  Kw, 


and  the  Klinkenberg ’s  constant  b . These  three  parameters  taken  in  conjunction  can 
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be  used  to  fully  characterize  the  single-phase  flow  of  either  liquid  or  gas  through 
concrete.  Consequently,  fluid  flow  rates  computed  by  using  these  three  intrinsic  flow 
parameters  of  concrete  are  significantly  different  from  those  predicted  by  using  a 
single-valued  permeability  in  a transient  heat  and  mass  transfer  analysis.  Pore 
pressure  development  (in  locations  and  time)  has  been  found  to  be  affected  and 
determined  mainly  by  the  assessment  of  the  fluid  flow  rates. 

The  flows  that  occur  in  a partially  saturated  concrete  exposed  to  nonisothermal 
conditions  involve  two  distinct  fluid  phases  (e.g.  liquid  water  and  gaseous  steam) 
permeating  simultaneously  through  the  porous  network.  Since  the  presence  of  the 
liquid  phase  can  reduce  the  permeability  of  the  concrete  to  the  gas  phase  and  vice 
versa,  we  employ  a postulate  commonly  used  to  allow  for  this  interaction  as 

krg=Rpg(n’ Sg)  (4.18) 

kn*  = RPM>  SJ 

where  krj!  stands  for  the  relative  permeability  (RP)  of  the  concrete  to  phase  n at  the 
degree  of  fluid  saturation  Sn , n represents  the  porosity,  and  the  functions  RP^(«,  S^) 
are  called  relative  permeability  functions  (RPF),  which  obey  the  bounds 
0<RP >,SJ<1. 

Using  the  relationship  S = 1 - Sw , a RP  of  concrete  to  gas  has  been  developed 

as  a function  of  the  degree  of  liquid  water  saturation  ( Sw ) and  an  empirical  parameter 

(A)  determined  by  the  porosity  n of  concrete.  Experimental  data  (Whiting  1988, 
Cabrera  and  Ghoddoussi  1 994,  Baroghel-Bouny  and  Chaussadent  1 996)  and  sorption 


isotherm  concept  (Hansen  1989)  were  used  by  the  author  to  develop  the  following 
relative  permeability  function 
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i?Pg(S„,/L)  = 10^-10^ 


(4.19) 


where 


A = 0.05- 22. 5« 


(4.20) 


Similarly,  for  a RPF  of  concrete  to  liquid  water,  the  ratio  of  the  water  permeability  to 
the  gas  permeability  has  been  combined  with  an  assumed  mirror  image  of  the  RPF  for 
gas  to  obtain 


Due  to  the  considerable  difficulties  involved  in  performing  and  interpreting  two- 
phase  (water-air)  RP,  very  little  data  is  available  to  construct  the  RP  of  concrete  to 
water.  However,  a remarkable  resemblance  in  the  RP-saturation  relations  in 
two-phase  (air-water)  and  three-phase  (air-water-oil)  systems  of  soils  is  found  in  both 
shape  and  variation  with  respect  to  water  saturation  (Lenhard  and  Parker 1987, 
Burdine  1953,  Parker  et  al.  1987).  Experimental  data  from  the  literature  (Burdine 
1953,  Parker  et  al.  1987)  also  indicate  that  the  proposed  trend  found  between  porosity 


RPW(SJ  - ctRPgi}  — Sw) 


(4.21) 


where 


(4.22) 
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and  sensitivity  of  RP  to  Sw  is  consistent  with  that  found  in  pore  size  distribution 
varying  with  concrete  types  (Figure  4-2). 
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Figure  4-2  Relative  permeability  functions  for  the  gas  phase 
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Concrete  as  a porous  medium  is  assumed  to  have  a similar  RP-saturation  relation  to 
that  of  soils.  Based  on  this  assumption,  the  shape  of  the  curve  of  RP  to  gas  can  be 
mirrored  to  construct  that  for  water.  However,  the  author  recognizes  limitation  of  the 
proposed  RPW  (1)  while  the  decrease  in  relative  air  permeability  is  more  pronounced, 
relative  water  permeability  may  mildly  increase  in  concrete  with  water  saturation,  and 
(2)  the  residual  degree  of  saturation  is  not  considered  in  the  proposed  RPW  whereas 

relative  water  permeability  of  soil,  experimentally  measured,  becomes  zero  if  Sw 


becomes  lower  than  the  residual  degree  of  saturation. 
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The  relative  gas  permeability  curves  generated  by  the  RPg  function  above  are 

plotted  for  various  values  of  porosity  in  Figure  4-2  along  with  the  data  points  used  in 
the  development.  Figure  4-3  presents  a comparison  between  Eqn.  (4.19)  and 
experimental  results  presented  by  Jacobs  (1998). 


By  conducting  gas  permeability  tests  on  concrete  specimens  at  various  states  of 
partial  liquid  saturation,  Jacobs  determined  RP  values  for  gas  flow  (i.e.  values  of  RPg 

as  a function  of  Sw ).  Upper  and  lower  bound  curves  for  Jacob’s  results  for  concrete 

specimens  that  had  W/C  ranging  approximately  from  0.45  to  0.8,  are  shown  in  Figure 
4-3.  For  comparison,  the  range  of  values  predicted  by  Eqn.  (4.19)  for  concrete  mixes 
having  W/(C+P)  ranging  approximately  from  0.3  to  0.6,  is  also  shown  as  a hatched 
region  in  Figure  4-3.  While  the  two  sets  of  curves  are  not  identical,  their  general  form 
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is  in  good  agreement.  Since  Eqn.  (4.19)  has  been  developed  using  experimental  data 
sources  other  than  Jacobs’  results,  the  agreement  shown  in  Figure  4-3  suggests  a 
degree  of  validity  in  the  relationships  presented  both  herein  and  by  Jacobs  (1998). 

4.4.2  Slip  flow  constant  of  concrete 

In  very  low  permeability  materials  such  as  HSC,  the  slip  flow  phenomenon 
known  as  Klinkenberg’s  effect,  can  be  dominant.  As  Eqn.  (4.15)  indicates, 
permeability  of  concrete  to  gas  flow  can  be  much  larger  than  the  intrinsic  gas 
permeability  k g if  a material  constant  b (i.e.,  Klinkenberg’s  constant)  is  significant. 

A functional  relationship  between  b and  the  gas  permeability  Kg  of  the  concrete  in 

Eqn.  (4.15)  as  shown  in  Figure  4-4  has  been  developed  based  on  (1)  permeability  test 
data  published  by  Whiting  (1988)  and  Klinkenberg  (1941)  and  (2)  based  on 
relationships  between  the  gas  permeability  and  the  water  permeability  for  concrete 
presented  by  Dhir  et  al.  (1989) 

£ _ (-0.5818  ln(A!g)  - 19.1213)  (4.23) 

where  Kg  is  in  units  of  m2  and  b is  in  units  of  atmospheres. 

Dhir’s  data  were  measured  from  cover  layers  of  concrete,  which  might  not  be 
capable  of  characterizing  bulk  concrete.  However,  the  primary  concern  in  this 
dissertation  is  on  pressure  and  temperature  profiles  near  the  surface  of  the  concrete 
where  thermal  spalling  occurs.  Thus,  Dhir’s  relationship  may  not  be  universally  used 
for  characterizing  concrete  at  all  depths  of  structural  elements,  but  it  is  well  suited  for 


the  situations  considered  here. 
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Figure  4-4  Slip  flow  constant  vs.  the  gas  permeability  of  concrete 


For  comparison,  the  slip  flow  values  predicted  by  Bamforth’s 
equation  (Bamforth  1987)  with  data  only  for  normal  weight  concrete  mixes  having 
W/(C+P)  ranging  approximately  from  0.33  to  0.41  and  Agg/C  ranging  approximately 
from  0.37  to  0.40  are  also  shown  in  Figure  4-4.  One  can  see  favorable  agreement  in 
the  range  of  a concrete  mixture  having  both  a low  W/C  and  a high  W/C.  For 
modeling  purposes,  an  approximate  range  of  slip  flow  constants  (shown  as  a hatched 
region  in  Figure  4-4)  is  proposed 


b - (1.635E-8)»(A’g.)-0'5227 


(4.24) 
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where  Eqn.  (4.24)  was  developed  by  Bamforth  (1987)  and  used  as  the  lower  bound  in 
Figure  4-4.  For  comparison,  an  equation  for  soils  given  by  the  American  Petroleum 
Institute  (1952)  is  also  plotted. 

4.4.3  Volume-averaged  quantities  for  thermal  properties  of  concrete 

Thermal  conductivity  and  diffusivity  of  concrete  are  generally  affected  by 
changes  in  the  type  and  quantity  of  aggregates.  The  concrete  considered  in  this  study 
is  assumed  to  be  made  of  siliceous  type  of  aggregates  such  as  quartzite  or  sand  stone 
with  Agg/(C+P)  = 3.7  to  6.2  and  W/(C+P)  - 0.26  to  0.6. 

Based  on  a spatial  volume  averaging  theory,  the  volume-averaged  thermal 
conductivity  used  in  Eqn.  (4.1 1),  which  corresponds  to  real  laboratory  measurable 
quantities,  are  assumed  to  be  a linear  function  of  the  degree  of  water  saturation  only 
(Pruess  1991) 


where  rc  (Sv)  represents  the  volume-averaged  thermal  conductivity.  Thermal 

conductivity  of  concrete  is  assumed  to  be  a relatively  weak  function  of  temperature, 
but  may  be  affected  significantly  by  moisture  content. 

However,  the  author  recognizes  that  temperature  significantly  influences  other 
properties  of  concrete  (Bazant  and  Kaplan  1996).  Particularly,  an  increase  in  porosity 
and  the  extent  of  dehydration  are  mainly  responsible  for  a sudden  jump  in  magnitude 
of  permeability  and  a noticeable  decrease  in  conductivity  of  concrete  at  elevated 
temperatures.  Bazant  (1997)  explained  that  at  elevated  temperatures,  the  micropores 
such  as  gel  pores  of  molecular  dimensions  that  are  inaccessible  in  an  ambient 


(4.25) 
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condition  become  accessible  to  the  fluids.  Rapid  increase  of  the  permeability  (by  an 
order  or  two  of  magnitude  increase)  can  be  observed  without  noticeable  change  of  the 
porosity.  Also,  microcracks  initiated  by  differential  thermal  expansion  occurring  at 
the  interfaces  between  cement  paste  (which  shrinks  as  it  is  heated)  and  aggregates 
expanding  during  fire  will  cause  opening  of  more  continuous  channels  for  fluid  flow, 
which  can  result  in  a substantial  increase  of  the  porosity.  These  phenomena  can  also 
contribute  to  a significant  increase  in  the  permeability  of  concrete  to  both  moist  air 
and  liquid  water  (Samaha  and  Hover  1992). 

More  importantly,  the  effects  of  moisture  movement  on  the  permeability  and 
conductivity  of  concrete  are  considered  in  modeling  simultaneous  flow  of  heat  and 
fluids.  For  example,  relative  permeability  to  moist  air  in  the  hydro-thermal  processes 
involved  in  the  formation  of  pressure  driven  moisture  layers  known  as  “ moisture  clog 
formation ” (Harmathy  1964),  is  coupled  with  thermal  conductivity  of  concrete 
between  quasi-saturated  and  de-saturated  zones. 

Thus,  the  presented  model  can  predict  the  fluid  phases’  involvement  with  heat 
transfer  that  includes  (1)  the  fluids  (moist  air  and  liquid  water)  being  driven  through 
concrete  by  pressure  gradient  and  being  accompanied  by  heat  convection  and  (2) 
evaporation  of  free  water  content  associated  with  latent  heat  of  water.  Hydro-thermal 
properties  used  here  for  structural  steel  and  two  types  of  concrete  are  obtained  from 
the  field  test  data  reported  in  literature  (Dhir  et  al.  1989,  Burg  and  Ost  1994,  Scanlon 
and  Mcdonald  1994)  and  are  given  in  Table  4-1.  Ongoing  experimental  research  in 
measuring  concrete  permeability  will  be  used  to  extend  Eqn.  (4.19)  to  include  the 
temperature  dependency  of  the  RP.  In  addition,  numerical  implementation  of  thermal 
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conductivity  of  concrete  as  a function  of  temperature  should  be  included  in  the  future 
research. 


Table  4-1  Material  properties  of  concrete  and  steel  used  in  numeric  simulations 


Material  Property 

Structural  Steel 

Concrete  I 

Concrete  11 

Porosity 

N/A 

13.3% 

8.2% 

Conductivity  (wet) 

52.0  W/m-C 

2.9  W/m-C 

2.6  W/m-C 

Conductivity  (dried) 

52.0  W/m-C 

1.7  W/m-C 

2.0  W/m-C 

Expansivity 

12.0E-6  1/C 

6.5E-6  1/C 

8.5E-6/1C 

Emissivity 

N/A 

0.8 

0.88 

Saturated  unit  weight 

7850  kg/  m3 

2200  kg/  m3 

2400  kg/  m3 

Specific  heat 

209  J/kg-C 

921  J/kg-C 

921  J/kg-C 

Intrinsic  gas  permeability 

N/A 

2.24E-17  m2 

6.34E-19  m2 

Intrinsic  liquid 
permeability 

N/A 

8.49E-19  m2 

4.82E-21  m2 

Klinkenberg's  Constant 

N/A 

2.44  MPa 
(24.1  atm.) 

19.45  MPa 
(192.0  atm.) 

28  Day  Compressive 
Strength 

N/A 

60  MPa 

90  MPa 

By  having  established  the  key  material  parameters  and  constitutive  equations 
used  in  the  governing  flow  equations,  several  example  simulations  will  be  presented 
to  study  thermodynamic  states  developed  in  concrete  exposed  to  fire.  In  the  following 
sections,  we  will  first  investigate  initial  degree  of  saturation  of  concrete  that  plays  a 
critical  role  in  prediction  of  pore  pressure  and  temperature  at  high  temperatures.  Then 
we  will  validate  the  numerical  model  developed  in  this  study  by  experimental  test 
results  reported  by  Consolazio  et  al.  (1998)  before  presenting  example  simulation 
results. 


4.5  Calculation  of  Initial  Conditions 

Initial  conditions  for  example  problems  presented  in  this  study  are  assumed  to 
be  uniform  conditions  of  pressure,  temperature,  and  saturation  throughout  the  entire 
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flow  domain  at  the  initial  state.  Initial  pressure  and  temperature  inside  the  flow 
domain  is  assumed  to  be  equal  to  the  atmospheric  pressure  and  a temperature  of  20 
°C,  respectively. 

Variability  of  environmental  conditions  affecting  concrete  in  situ  and  a variety 
of  concrete  mixtures  are  taken  into  consideration  to  determine  the  initial  conditions  of 
concrete  at  ambient  conditions.  Initial  degrees  of  saturation  were  evaluated  for  two 
types  of  concrete  (Table  4-1)  concrete  I with  W/C=0.5-0.6  and  Agg/C-5. 0-6.0  and 
concrete  II  with  W/(C+P)=0.28-0.3  and  Agg/(C+P)=3.8-4.2,  in  an  ambient  condition 
of  60  % RH  at  a temperature  of  20  °C. 

The  moisture  content  in  a partially  saturated  concrete  can  be  computed  using 
sorption  isotherms  of  concrete  as  a percentage  of  the  free  moisture  content.  Using  the 
porosity,  W/(C+P),  Agg/(C+P),  and  a saturated  unit  weight  of  the  concrete,  the  water 
content  of  the  saturated  concrete  is  computed  and  subsequently  converted  into  a 
“water  content  by  mass  of  dried  hardened  cement  paste  (HCP)”  (Baroghel-Bouny  and 
Chaussadent  1996).  This  step  of  converting  the  “normal”  mass  water  content  into  one 
that  is  referenced  with  respect  to  only  the  portion  of  the  concrete  mass  attributable  to 
the  cement  paste  was  performed  so  that  the  concrete  isotherms  of  Baroghel-Bouny 
and  Chaussadent  (1996)  could  be  used. 

Taking  saturated  concrete  having  Sw  - 1 .0  and  oven-dried  (until  less  than  3% 
change  in  weight  is  observed)  concrete  having  Sw  = 0.0 , an  intermediate  value  of 

degree  of  saturation  can  be  computed  as  a percentage  of  the  saturated  moisture 
content  by  using  the  weight  differences  between  the  saturated  and  dry  concrete 
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Sii=_^*i*se.xioo  <4'26) 

wsat;hcp 

where  Sw  represents  the  degree  of  saturation,  wsarhcp  = water  content  (by  mass  of 

HCP)  of  concrete  in  a saturated  condition,  and  wrh.hcp  = water  content  when  the 

concrete  is  air-dried  at  a particular  RH  level  at  the  temperature  of  20  °C. 

For  example,  consider  a concrete  sample  made  of  concrete  I in  Table  4-1  having 
a porosity  of  0.133  (13.3%),  a W/(C+P)=0.5,  an  Agg/(C+P)=5.8,  and  a fully  saturated 
water  content  of  wsat:hcp=  0.1989  (19.89%).  Assume  that  the  concrete  member  has 

reached  hydro-thermal  equilibrium  at  approximately  60-65%  relative  humidity  at  the 
temperature  of  20  °C.  Isotherm  data  from  experiments  by  Baroghel-Bouny  and 
Chaussadent  (1996)  can  be  used  to  find  that  the  water  content  is  approximately  0.334 
(33.4%)  by  mass  of  HCP.  Using  Eqn.  (4.26)  and  the  concrete  isotherms,  a degree  of 
initial  saturation  for  the  NSC  is  approximated  as  Sw  = 0.4  in  a range  of  60-65%  RH. 
Similarly,  for  a concrete  sample  having  a porosity  of  0.082  (8.2%),  a W/(C+P)=0.28, 
an  Agg/(C+P)=3.8,  and  a fully  saturated  water  content  of  wsaVhcp  = 0.10  (10%),  a 

degree  of  initial  saturation  for  the  higher  porosity  concrete  is  approximated  as  Sw  = 
0.6  in  a range  of  60-65%  RH. 

4.6  Comparison  of  Numerical  Results  with  Experimental  Data 

Before  conducting  simulations  by  the  proposed  modeling  techniques,  validation 
against  published  data  for  saturated,  but  unreinforced,  mortar  exposed  to  severe 
heating  are  performed.  A 2-D  simulation  is  performed  and  results  are  compared  to  the 
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experimental  data  reported  by  Consolazio  et  al.  (1998).  The  original  2-D  TOUGH 
mesh  used  by  Consolazio  et  al.  (1998)  in  numerical  simulations  was  acquired  and 
lightly  modified  for  TOUGH2  simulations  presented  here.  Three  new  features  are 
included,  (1)  an  “ effective ” conductive  boundary  approach,  a modeling  technique 
used  for  a radiation  boundary  condition  between  a specimen  and  a radiant  heat 
source,  (2)  the  proposed  relative  permeability  function , and  (3)  the  relationship 
between  the  permeability  and  slip  flow  constant. 

Because  water  permeability  of  the  cement  mortar  specimen  considered  by 
Consolazio  et  al.  was  not  available,  a logarithmic  relationship  between  the  intrinsic 
gas  permeability  data  and  the  intrinsic  water  permeability  data  of  concrete  from 
experimental  data  of  Dhir  et  al.  (1989)  is  used  to  approximate  the  water  permeability 
of  a cement  mortar  whose  intrinsic  gas  permeability  is  8.3424E-16  m 2 


Kg  = exp  (-9.626  + 0.69  ln(£J)  { j 

However,  concern  over  possible  deviation  of  the  approximated  value  from  the  actual 
water  permeability  due  to  the  type  of  aggregates  used  and  the  Agg/(C+P)  were 
believed  to  significantly  affect  the  water  permeability.  Therefore,  upper  and  lower 
bounds  of  the  ratio  a (i.e.,  the  ratio  of  the  intrinsic  water  permeability  to  the  intrinsic 
gas  permeability  of  a concrete)  are  used  for  simulations  as 


0.038  <a  < 0.068 


(4.28) 


The  upper  bound  of  the  ratio  a was  computed  by  Eqn.  (4.27)  whereas  the  lower 
bound  value  was  extrapolated  for  a concrete  having  a 14.5  % porosity  from  the  data 
available  for  concrete  having  porosity  in  a range  between  7.5%  and  13.3%  in 
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Whiting’s  experiments  (1988).  The  results  of  simulating  the  conditions  considered  by 
Consolazio  et  al.  (1998)  are  shown  in  Figure  4-5. 

As  illustrated  in  Figure  4-5a,  pore  pressures  occurring  at  a position 
approximately  1 9 mm  from  the  heated  surface  of  the  model  are  compared  to 
experimental  data  measured  at  the  same  location.  The  pore  pressure  dip  for  the 
experimental  case  was  reported  as  attributable  to  voids  created  by  the  placement  of 
pore  pressure  transducers  in  testing,  “the  dip  is  a product  only  of  the  specific 
instrument  setup  and  not  a phenomenon  that  would  normally  occur  in  a non- 
instrumented  specimen"  (Consolazio  et  al.  1998),  and  thus  is  not  of  great  concern  in 
terms  of  validating  the  numeric  model.  More  importantly,  the  initial  peaks  for  the 
simulated  and  experimental  curves  match  in  magnitude  and  location  (time). 

One  may  note  that  the  peak  pore  pressures  predicted  by  the  numerical 
simulation  with  the  lower  bound  ratio  very  closely  match  the  actual  pore  pressure 
measured  experimentally.  Despite  the  same  volume  fraction  of  voids  (porosity),  a 
cement  mortar  specimen  can  have  a denser  pore-size  distribution  than  concrete.  This 
is  believed  to  result  in  lowering  the  water  permeability  of  the  cement  mortar 
compared  to  a concrete  having  the  same  porosity.  Thus,  the  model  with  the  lower 
bound  value  of  a seems  to  be  better  matched  with  pore  pressure  data  measured  by 
Consolazio  et  al.  (1998). 

Temperature  data  at  19  mm  depth  as  predicted  by  the  simulation  (Figure  4-5b) 
shows  an  overall  match  with  the  experimental  curve  although  the  simulated 
temperature  curve  at  1.5  mm  depth  diverges  significantly  from  the  experimental  data 


Temperature  (C)  ' — ' Pressure  (MPa) 
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Pore  pressures  measured  experimentally  and  predicted  using  a 2-D  simulation 

model 


Time  (sec) 


(b)  Temperatures  measured  experimentally  and  predicted  using  a 2-D  simulation 

model 


Figure  4-5  Comparison  of  numerical  prediction  to  experimental  measurement 
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This  indicates  that  conductive  and  convective  heat  transfer  through  the  specimen  is 
well  predicted  yet  the  temperature  changes  at  locations  nearer  to  the  surface  exposed 
to  the  heat  source  may  have  been  involved  with  errors  in  modeling  the  radiant 
boundary  condition  and  in  the  assumed  location  of  the  thermocouples  as  reported  in 
the  literature. 


4.7  Finite  Difference  Simulation 

Having  validated  the  modeling  proposed  techniques,  focus  will  now  be  turned 
to  studying  the  effects  of  reinforcing  steel  on  moisture  movement  in  heated  R/C 
columns.  2-D  simulation  models  are  constructed  for  a longitudinally  R/C  column 
made  of  concrete  I (Table  4-1)  with  a 0.38  x 0.38  m square  cross  section  exposed  to  a 
radiant  fire  as  shown  in  Figure  4-6a.  Due  to  symmetry  of  the  geometry  and  the 
thermal  loading  conditions,  partial  symmetry  models  were  used,  which  are  shown  in 
Figure  4-6b  and  4-6c  for  each  heating  condition  simulated. 

4.7.1  Modeling  boundary  conditions 

One-half  and  one-quarter  of  the  cross  section  were  discretized  by  symmetry. 
Along  the  symmetric  boundaries,  zero  flux  (no  flux  of  mass  or  heat)  boundary 
conditions  were  specified  whereas  the  surfaces  exposed  to  the  heat  source  were 
modeled  as  a radiant  heat  flux  boundary  condition.  For  the  radiant  heat  flux  boundary 
denoted  by  T? , heat  source  is  modeled  using  ASTM  El  19 

= on  r,  (4.29) 

where  Kebound  represents  the  effective  conductivity  on  the  boundary  between  the 
radiant  heat  source  and  the  exposed  surface.  Particularly  for  the  case  of  one-side 
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thermal  loading  condition  (Figure  4-6  a),  the  boundaries  on  ambient  air  are  modeled 
as  a Dirichlet  boundary,  i.e.  a constant  pressure  and  a constant  temperature  condition, 
by  using  elements  that  are  very  permeable  and  have  very  large  volume  relative  to  the 
surface  elements  so  as  to  ensure  the  constant  thermodynamic  state  in  ambient  air 

T (x,  y, t)  = T0  (x, y ) = Tamb  on  T T (4 .30) 

P(x,y,t)  = P0(x,y)  = Palm  on 

where  Tamh  represents  ambient  temperatures  at  the  initial  state  on  the  temperature 
boundary  Tr  and  Palm  is  the  atmospheric  pressure  that  is  constant  on  the  pressure 
boundary  r;,  throughout  the  time  domain. 

4.7.2  2-D  simulation  results 

From  a 2-D  simulation  with  one-side  heating  boundary  condition,  variation  in 
build-up  of  pore  pressure  along  the  heat  flux  direction  is  observed  along  with 
discontinuity  in  development  of  thermal  gradients  due  to  presence  of  the  longitudinal 
bars  (Figure  4-7). 

A time  history  of  the  pore  pressure  build-up  is  plotted  in  Figure  4-7a,  which  is 
sampled  at  a location  of  x = 140  mm  (Figure  4-6b).  For  a better  resolution,  it  is 
plotted  for  only  the  first  120  mm  along  the  y-coordinate.  A quasi-saturated  zone  is 
developed  in  front  of  the  steel  bars  much  earlier  than  that  occurring  at  the  adjacent 
concrete  domain  at  x = 140  mm  (Figure  4-7b).  This  early  moisture  clog  formation 
affects  the  distributions  of  pore  pressure  and  temperature.  Most  notably,  the  peak  pore 
pressure  in  the  front  region  of  the  steel  bar  ( y = 38  mm ) reaches  up  to  4.35  MPa  after 
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(a)  2-D  geometric  dimensions 
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T(x,y,t)  = Tamb 
P(x,y,t)  = Pam 


Radiant  Heat  Flux 
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(b)  A 2-D  discretized  mesh  with  an  one-side  (c)  A 2-D  discretized  mesh  with  a four-side 
heating  condition  heating  condition 


Figure  4-6  2-D  simulation  models 


Heat  Flux  Boundary 
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(a)  Pore  pressure 


(b)  Degree  of  saturation 


(c)  Temperature 


Figure  4-7  Thermodynamic  state  variables  determined  by  a 2-D  simulation 


after  30  minutes  of  heating  and  the  peak  of  pore  pressure  at  a point 
( x = 1 00 , = 38  mm ) is  predicted  as  less  than  3 . 5 MPa. 

A discontinuity  in  the  temperature  field  occurs  at  the  front  region  of  the  steel 
reinforcing  bar  (Figure  4-7c).  Because  of  the  early  formed  moisture  clog  caused  by 
the  impermeable  steel,  heat  carried  by  the  moisture  flow  in  the  surface  concrete  is 
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accumulated  whereas  heat  can  transfer  quickly  through  the  steel  bar,  thermal 
diffusivity  of  steel  is  greater  than  that  of  concrete.  Therefore,  a steeper  thermal 
gradient  is  observed  in  the  concrete  region  in  front  of  the  steel  whereas  almost  no 
thermal  gradient  is  developed  within  the  steel  layer.  Consequently,  the  slope  of 
temperature  field  (thermal  gradients)  noticeably  changes  as  a quasi-saturated  zone 
forms  and  affects  the  thermally  driven  moisture  flow. 

Figure  4-7b  also  indicates  that  a quasi-saturation  zone  is  formed  behind  the  steel 
reinforcement  whereas  the  degree  of  saturation  at  the  front  decreases.  The  build-up  of 
pore  pressure  in  the  front  of  the  steel  becomes  large  enough  to  force  the  clogged 
moisture  to  permeate  through  the  concrete  area  around  the  steel.  Also,  the  degree  of 
saturation  in  front  of  the  steel  decreases  as  evaporation  of  the  water  content  in  the 
clogged  region  begins.  The  more  temperature  and  pressure  fields  rise,  the  more 
evaporation  occurs  in  the  region.  Eventually  the  pressure  build-up  dissipates  when  the 
region  is  desaturated.  Apparent  permeability  of  concrete  to  gas  is  thought  to  be  the 
reason  for  quick  dissipation  because  the  intrinsic  permeability  to  moist  air  is  higher 
than  that  to  liquid  water  due  to  the  slip  flow  phenomenon.  Thus,  peak  of  pore 
pressure  build-up  is  found  to  be  in  the  moisture-clogged  regions  rather  than  the 
desaturated  zones.  Therefore,  presence  of  the  steel  reinforcement  impedes  moisture 
movement  through  the  concrete  area.  This  clogging  results  accumulation  of  moisture 
in  front  of  the  steel. 

4.7.3  Comparison  of  2-D  simulations  modeled  by  various  fire  conditions 

In  this  section,  three  different  fire  conditions  are  simulated.  Although  it  is 
known  that  a higher  rate  of  heating  rapidly  induces  a steeper  thermal  gradient  and 
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larger  pore-pressure  build-up  compared  to  a slow  rate  of  heating,  qualitative 
computation  of  pore  pressure  and  its  time  history  were  of  interest  here.  For  example, 
the  relationship  between  the  intensity  of  the  heating  rate  and  time  histories  of 
moisture-clog  formation  and  temperature  distribution  in  the  presence  of  steel 
reinforcement  need  to  be  adequately  addressed  for  a better  understanding  fire 
resistance  of  high-strength  concrete. 

A 2-D  simulation  of  all  four  sides  of  the  reinforced  column  exposed  to  fire  is 
performed.  The  data  of  pore  pressure  and  saturation  are  plotted  in  Figure  4-8  along 
the  line  of  symmetry  as  indicated  in  Figure  4-6c.  Comparing  to  a 2-D  simulation  with 
one-side  heating  boundary  condition,  earlier  formation  of  moisture  clogging  in  the 
four  side  heating  simulation  is  observed  at  t=l  500  sec.  (25  minutes)  and  pore  pressure 
at  x = 100  mm  and  y = 38  mm  is  predicted  as  6.0  MPa.  The  more  heat  comes  into 
the  system,  the  more  rapid  formation  of  moisture  clog  is  induced  and  the  greater  pore 
pressure  builds  up.  Therefore,  the  influence  of  moisture  movement  in  pore  pressure 
build-up  is  more  noticeable  in  the  simulation  of  four  side  heating  condition. 

Next,  a 2-D  simulation  for  the  case  of  all  four  sides  of  the  column  exposed  to  a 
hydrocarbon  pool  fire  (modeled  by  ASTM  El  529)  is  performed.  At  a geometric 
location  (y'  = 38  mm , see  Figure  4-9),  the  pore  pressure,  saturation,  and  temperature 
data  are  plotted  and  compared  with  other  simulation  results  predicted  under  different 
fire  conditions  as  shown  in  Figure  4-9.  The  case  of  ASTM  El  529  heating  shows  the 
most  pronounced  influence  of  movement  of  the  saturation  front  in  pore  pressure 
because  the  hydrocarbon  fire  has  a much  steeper  temperature  rise  in  the  first  ten 


minutes. 
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Therefore,  the  rate  of  heating  and  boundary  heating  condition  significantly 
affects  development  of  thermodynamic  states  of  concrete,  the  pressure  and 
temperature  profiles,  particularly  during  the  first  twenty  minutes  of  simulation  time. 
The  presence  of  steel  reinforcement  produces  more  influence  in  pressure  and 
temperature  fields  in  concrete  when  exposed  to  high  rate  of  heating  than  when 
exposed  to  low  rate  of  heating. 

4.7.4  3-D  simulation  results 

Effects  of  steel  reinforcement  on  heat  and  mass  transfer  in  the  concrete  domain 
have  been  found  to  be  severe,  particularly  in  high  rate  heating.  Even  a larger  build-up 
of  pore  pressure  and  a significant  change  in  temperature  distribution  (development  of 
steeper  temperature  gradients)  are  predicted  by  simulating  reinforced  steel  in  concrete 
structural  elements  exposed  to  fire. 

A 2-D  simulation  model  including  shear  reinforcement  would  not  be 
appropriate  simply  because  shear  reinforcement  would  completely  prevent  the 
moisture  flow  from  reaching  a deeper  region  beyond  the  impermeable  layer  of  steel 
(Figure  4- 10a  and  Figure  4- 10b).  Instead,  a 3-D  simulation  model  is  constructed  that 
is  capable  of  simulating  three-dimensional  fluid  flow  through  concrete  surrounding 
reinforcing  steel  bars  including  both  longitudinal  and  shear  reinforcement. 

Due  to  symmetry  of  the  geometry  and  the  thermal  loading  condition  (modeled 
using  ASTM  El  529),  only  one-eighth  of  the  actual  column  segment  is  discretized 
with  no  flux  boundary  conditions  at  all  the  symmetric  surfaces  (Figure  4- 10c).  A 
large  rise  of  pore  pressure  is  observed  (Figure  4-1  la)  at  a depth  of  z’=  5 mm  along 
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9.2  MPa 


Figure  4-8  Variation  of  pore  pressures  and  degrees  of  saturation 
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Simulation  with  one  side  heating  by  El  19  

■ Simulation  with  four  side  heating  by  El  19  

/ 

/ V / 

Simulation  with  four  side  heating  by  El 529  

/ 

f I i 

/ 

✓ 

/ 

;✓ 

_ ✓ ’ 

0 200  400  600  800  1000 

Time  (sec) 


Figure  4-9  Comparison  of  2-D  simulation  results 
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Longitudinal  bars 


(a)  A schematic  sketch  of  a reinforced  concrete  column  segment 


Radiant  Heat  Flux 


r 

t \ 

- »■> 

L _N 

f 

L \ 

Shear  reinforcement 


9 

/ 

Longitudinal 

reinforcement 


f 


380  mm 


I I I 


-I 


"K> 

4^—10 


4^ 

O 

4= — 10 

to 
-o 

3 
3 


(b)  Geometric  dimensions 
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(c)  A discretized  mesh 


Figure  4-10  Geometry  of  a 3-D  simulation  mesh  and  a thermal  loading  boundary 


Ill 


Figure  4-1 1 Variations  of  pore  pressures  and  degrees  of  saturation 


the  y’-axis.  A quasi-saturated  zone  (Figure  4-1  lb)  is  rapidly  developed  in  concrete 
regions,  occurring  at  locations  adjacent  to  / = 20  mm  after  5 minutes  of  heating. 

This  rapid  formation  of  a moisture  clog  induces  the  peak  pore  pressure  ( y'  = 30  mm ) 
that  exceeds  8 MPa  after  7 minutes  of  heating.  In  particular,  a larger  pressure  build-up 
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(5.5  MPa)  was  predicted  in  the  3-D  simulation  along  the  boundaries  between  the 
concrete  regions  and  the  shear  steel  at  Point  A-l  in  Figure  4-1 1 whereas  pressure  less 
than  1 .5  MPa  was  predicted  ( y'  = 60  mm , z'  = 15  mm ) in  the  2-D  simulation  after  10 
minutes.  The  simulation  results  indicate  that  impermeable  steel  reinforcing  bars  alter 
moisture  flow  in  the  concrete  matrix.  Consequently,  the  temperature  distribution 
changes  and  very  high  pore-pressure  develops  in  the  concrete. 

Two  important  causes  for  explosive  thermal  spalling  of  concrete,  i.e.,  thermal 
gradients  and  pore  pressures,  are  found  to  be  significantly  altered  by  the  presence  of 
steel  reinforcement  in  heated  high-strength  concrete  structural  members.  From  the 
simulation  results,  very  high  pressure  build-up  can  be  developed  in  reinforced 
concrete  structures  made  of  very  low  permeable  concrete  such  as  concrete  with  a 28 
day  compressive  strength  of  more  than  60  MPa  with  a low  W/C  ratio. 

In  the  following  chapter,  the  discussion  will  be  expanded  to  include  predictions 
not  only  of  pore  pressure  and  temperature,  but  more  importantly  the  stresses  produced 
by  these  effects.  The  effects  of  steel  reinforcement  in  the  coupled  moisture  and  heat 
transfer  phenomena  will  be  addressed  in  Chapter  6. 


CHAPTER  5 

FINITE  ELEMENT  STRESS  ANAYSIS  OF  A REINFORCED  HIGH-STRENGTH 
CONCRETE  COLUMN  IN  SEVERE  FIRES 


The  objective  of  this  chapter  is  to  develop  finite  element  stress  analysis  models, 
which  can  predict  stress  states  developed  in  a reinforced  concrete  structural  element 
exposed  to  rapid  heating.  Obtained  from  finite  difference  models  of  simulating  coupled 
heat  and  mass  transport  phenomena  in  heated  reinforce  concrete  elements,  transient 
temperature  profiles  are  used  as  prescribed  boundary  conditions  for  subsequent  finite 
element  thermoelastic  stress  analysis.  Thus,  a computational  methodology  is  presented  to 
analyze  stress  states  associated  with  coupled  hydro-thermal  behavior  of  concrete  exposed 
to  severe  fire  conditions.  Multidimensional  finite  element  models  that  characterize 
structural  response  of  reinforced  concrete  columns  exposed  to  radiant  heating  conditions 
are  developed.  The  principle  of  superposition  and  a volume-averaging  theory  are  used  to 
compute  thermally  induced  effective  stresses  that  are  potentially  associated  with  thermal 
spalling  of  high-strength  concrete.  Effects  of  steel  reinforcements  on  coupled  heat  and 
moisture  transport  phenomena  are  also  considered. 

5.1  Goals  and  Scope 

In  the  past,  explosive  spalling  of  concrete  structures  has  been  investigated  by 
making  use  of  numerical  models  that  analyze  mass  and  heat  transfer  phenomena  in  plain 
concrete  exposed  to  fire.  Kodur  (1999)  reported  that  thermal  spalling  of  intensively 
heated  concrete  would  result  from  pore  pressure  and  thermal  stress  development  in 
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concrete  structural  elements.  In  particular,  high  strength  concrete  (HSC)  exposed  to 
rapid  heating  has  a greater  tendency  to  explosively  spalling  than  does  normal  strength 
concrete  (NSC)  (Phan  and  Carino  1998,  Sullivan  2001).  Previous  studies  by  Kodur 
(1997)  and  Kalifa  et  al.  (2000)  have  theorized  that  the  lower  permeability  of  HSC — 
compared  to  that  of  NSC — is  responsible  for  the  higher  susceptibility  of  HSC  to 
violent  spalling. 

In  past  published  articles  (Ymazaki  et  al.  1995,  Moyne  et  al.  2000),  no  major 
efforts  in  modeling  of  a coupling  between  thermo-mechanical  stress  and  pore 
pressure  have  been  found.  Only  a limited  number  of  computational  analyses  that 
focused  on  hydro-thermal  behavior  of  concrete  at  high  temperatures  have  been 
recently  reported  (Gawin  et  al.  1999,  Tenchev  et  al.  2001).  Consequently,  very  few 
simulation  codes,  which  are  capable  of  simulating  the  coupled  hydro-thermal  and 
structural  behavior  of  concrete  structures,  have  been  developed  in  this  research  field. 
The  needs  for  the  development  of  a new  numerical  tool  capable  of  simulating  such 
coupled  processes  are  thus  evident  (Phan  and  Carino  2000).  The  numeric  tool 
developed  here  will  assist  researchers  in  gaining  insight  into  the  thermo-mechanical 
behavior  of  concrete  structures  during  a transient  process  of  coupled  heat  and  mass 
transfer.  Implementation  of  the  new  analytical  model  is  accomplished  by  adding  new 
thermal  modeling  features  to  the  existing  PLASFEM  nonlinear  finite  element  code 
developed  by  Dr.  Michael  McVay  of  the  University  of  Florida,  Department  of  Civil 
and  Coastal  Engineering,  Geotechnical  Engineering  group. 

This  chapter  also  focuses  on  qualitatively  understanding  the  development  of 
stress  states  in  concrete  exposed  to  coupled  heating  and  moisture  migration  processes. 
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To  conduct  the  study,  a three-dimensional  combined  finite  difference  (TOUGH2)  and 
finite  element  model  (PLASFEM)  has  been  developed  that  is  capable  of  analyzing  the 
transient  heat  conduction  in  partially  saturated  concrete  and  predicting  states  of 
effective  stress.  An  existing  commercial  program,  ADINA  (ADINA  2002)  is 
employed  to  model  a 2-D  axisymmetric  thermoelastic  problem  and  to  validate  the 
present  3-D  finite  element  model.  The  coupling  effects  of  moisture  flow  on  the 
development  of  the  thermodynamic  temperature  field  are  studied  in  consideration  of 
the  following  important  characteristics  of  heat  and  mass  transfer  through  concrete 

1 . Strong  coupling  of  pore  pressure,  temperature,  and  moisture  content. 

2.  Multi-phase  mass  transfer  (liquid  and  gaseous  phases). 

3.  Nonlinear  characteristics  of  key  material  parameters  such  as  permeability, 
conductivity,  density,  and  porosity  and  their  dependency  on  pressure,  and 
degree  of  saturation. 

5.2  Governing  Differential  Equations  for  the  Mass  and  Energy  Balance 

The  balance  laws  used  in  this  study  treat  concrete  as  a multiphase  continuum 
analogous  to  that  of  a simple  continuum.  Global  balances  apply  to  the  multiphase 
continuum  as  a whole  via  an  averaging  process,  i.e.,  volume  averaging,  which  is  used 
for  a framework  to  construct  mathematical  formulations  of  governing  field  equations 
at  the  macroscopic  level,  the  first,  the  mass  balance  law  is  reviewed  for  a multi-  phase 
mixture,  which  is  a collection  of  overlapping  continua  called  phase. 


— + div pv  - 0 
8t 


(5.1) 
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where  each  phase  n occupies  a fraction  of  the  overall  mixture  volume  defined  by  a 
volume  fraction  Tj*(x,t)  and  has  the  volume-averaged  mass  density  p with  the  mean 
velocity  v(x,/)  defined  in  Chapter  4.  The  first  term  on  the  left  side  accounts  for  the 
accumulation  of  phases  and  is  known  as  the  mass  concentration  in  the  fluid  phases 
(water  and  water  vapor)  and  the  second  term  models  the  advection  and  diffusion  of 
the  fluid  phases  with  the  mean  macroscopic  flow  velocity. 

For  the  energy  balance,  we  have  a governing  equation  describing  the 
macroscopic  energy  balance 


where  the  terms  have  the  following  physical  interpretations  p — is  the  rate  of 

dt 

change  of  internal  energy,  V»q  is  the  rate  of  heat  flow,  -cr : Vv  represents  heating 
caused  by  compression  and  dissipative  momentum  transfer  (e.g.,  elastic  and  plastic 
deformation),  and  - ph  is  the  rate  of  heating  by  external  supplies. 


Based  on  the  Helmholtz  free  energy  and  the  Fourier  law  of  heat  conduction,  a 
coupled  heat  conduction  equation  that  describes  the  energy  balance  of  Eqn.  (5.2)  in  a 
concrete  member  can  be  expressed  at  the  macroscopic  level  (Hsu  1986) 


p — + V*q  - cr : Vv-p/?  = 0 
dt 


(5.2) 


dE 


5.3  The  Coupled  Heat  Conduction  Equation 


Pcp  — = (kTi)i+Q+T 
dt 


(5.3) 
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where  Cp  represents  a volume-averaged  specific  heat  at  a constant  pressure,  the 

absolute  temperature  of  a representative  element  volume  (REV)  (Lewis  and  Schrefler 
2000)  of  concrete  is  denoted  by  T , k represents  a volume-averaged  thermal 
conductivity  that  is  assumed  to  be  a function  of  temperature  and  degree  of  liquid 


water  saturation,  Q represents  heat  source  term  in  the  REV,  T ^ represents  a 


volume-averaged  thermo-elastic  coupling  factor.  In  this  conceptual  formulation  of 
thermal-energy  associated  with  mechanical  strain  energy,  only  a strictly  reversible 
energy  balance  is  considered  in  this  study. 

In  regard  specifically  to  the  mass  and  energy  balance  framework  equations 
described  above,  the  following  assumptions  are  made  1 ) the  concrete  material  can  be 
characterized  as  a continuum  via  a volume  averaging  process  such  that  all  physical 
quantities  (e.g.,  mechanical  deformation,  temperature)  must  follow  continuous 
functions,  2)  changes  in  material  temperature  due  to  mechanical  strain  are  negligibly 
small  and  can  therefore  be  neglected  (i.e.,  within  the  scope  of  conducting  the  heat  and 


moisture  flow  analyses  in  this  study,  the  quantity 


da 

fji  IJ 


en 
dT  ,J 


in  Eqn  (5.3)  is 


neglected). 


5.4  The  Principle  of  Stress  Superposition 

To  analyze  the  combined  effects  of  moisture  movement  and  temperature 
gradient  development  in  the  computation  of  thermally  induced  stresses,  the  principle 
of  superposition  is  used  in  solving  the  linear  momentum  balance  equation  with  mixed 
boundary  conditions.  Discussion  regarding  stress  states  in  heated  concrete  shall  be 
made  at  a macroscopic  level  where  continuous  distributions  of  the  phases  are 
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assumed.  This  means  that,  at  every  particle,  all  phases  are  present  at  the  same  time  in 
proportion  to  their  volume  fractions. 

To  integrate  microscopic  level  stresses  into  the  macroscopic  level  of  interest  of 
continuum  mechanics,  we  apply  the  volume-averaging  process  to  compute  stress 
components  of  a total  stress  state  at  a thermodynamic  equilibrium  state.  In  the 
following,  total  stress  at  the  macroscopic  level  is  defined  as  the  sum  of  averaged 
quantities  of  stress  components  through  summation  of  all  the  momentum  balance 
equations  written  for  the  solid  and  fluid  phases 


where  continuity  of  the  stress  field  across  the  solid-fluid  interface  is  assumed  to  be 
valid  so  that  the  linear  momentum  balance  equation  for  the  multi-phase  concrete 
medium  is  obtained. 

Therefore,  averaging  the  microscopic  stress  quantity  over  the  REV,  we  have 


where  q is  the  total  stress  tensor  macroscopically  averaged  over  the  REV  dVREV  . 
The  total  stress  tensor  consists  of  two  distinct  components 


V*  (rj'cr5  + 71*0-*)  + (77  V + 77*  p*  )g  = 0 


(5.4) 


(5.5) 


(5.6) 


where  dVs  and  dVf  represents  the  volumes  of  the  solid  phase  and  the  fluid  phase, 
respectively.  We  can  further  decompose  the  stress  tensor  of  the  fluid  phases  into  two 
parts  such  that 
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q dVREV  i ,yS  q dVm  + hr'  ^ dVf  hvw  q dVm  + dVf  ^yg  qdVm">  dV" 


(5.7) 


where  dV"  and  dVg  are  the  volume  of  liquid  water  and  the  volume  of  the  mixture  of 
air  and  water  vapor,  respectively. 


dVs  dV}  ,dVw  dVg  . 
- dVmv  ~s  dVBEU  dVf  ~v  dVf  ~g) 


(5.8) 


REV 


Using  the  volume  fraction  term  r)n  ( n - s,  w,  g solid,  liquid  water,  and  gas 
respectively),  Eqn.(5.8)  can  be  expressed  in  terms  of  the  porosity  n and  the  degree  of 
saturation  Sn  to  a particular  fluid  phase  n in  a partially  saturated  concrete 


q = {\-ri)qs+n-{Swqw  + Sgqg) 

where  qn  represents  the  stress  tensor  in  phase  n of  concrete.  If  shear  stresses  in  fluid 

phases  (liquid  water  and  the  air — water  vapor  mixture)  are  assumed  to  be  negligible, 
then  the  stress  tensors  of  the  fluid  phases  may  be  written  in  terms  of  hydrodynamic 
pressure  pK . We  can  write 


q = (l-n)qs-n- (Swpwl  + Sgpg\)  ^ 

Assuming  that  local  thermodynamic  equilibrium  exists,  then  pw  = pg  = p where  p 

is  the  internal  pore  pressure  (the  same  quantity  obtained  from  TOUGH2  analyses), 
and  the  total  stress  tensor  is  may  be  expressed  as 


q - (1  - n)  ■ qs  - n ■ pi 


(5.11) 
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Thus,  the  effective  stress  is  written 


=<Z  + n-Pl  = <I  + Pavgl  1 

where  pmg I represents  the  volume-averaged  hydrostatic  stress  tensor  at  a 
thermodynamic  equilibrium  state.  The  volume-averaged  total  stress  tensor  a then 
consists  of  mechanically  induced  stress  and  dilatation-induced  thermal  stress 


q"s  = q + n P I 


(5.13) 


where  q”  is  the  effective  stress  acting  on  the  solid  skeleton.  The  quantity  a" 
represents  the  stress  carried  solely  by  the  solid  skeleton  and  therefore  may  be  used  in 
assessing  the  potential  for  fracture  related  material  failure.  However,  in  evaluating 
effective  stress,  alternative  porosity  parameters  may  be  used  to  more  accurately 
reflect  the  failure  surface  involved.  Bazant  and  Kaplan  (1996)  note  that  boundary 
porosity,  np , may  be  more  appropriate  in  computing  effective  stresses  developed  in 


the  concrete  solid  skeleton.  They  define  np  as  the  porosity  of  the  “macroscopically 

planar  but  microscopically  tortuous  weakest  section  through  the  materials”  and  note 
that  np  = 0.9  is  a reasonably  representative  value  for  concrete.  Thus,  in  the  present 

study,  the  linear  momentum  balance  of  the  concrete  system  is  expressed  in  term  of 
the  effective  stress 


v.g;=o 


(5.14) 


Using  a generalized  Hooke’s  law  (Malvern  1969),  the  total  stress  tensor  can  be 
expressed 
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°u=c>jkt{ekt-Xkt(£n} 

= (K-ZgM,  + 2Ge'j  - 3Ka(AT)Su 


(5.15) 


where  Xkt  represents  a volume-expansion  coefficient  tensor,  K is  the  bulk  modulus, 

G is  the  shear  modulus,  is  the  volumetric  strain,  8 is  the  Kronecker  delta,  s'n  are 

the  deviatoric  strain  components,  a represents  a linear  expansion  coefficient,  and 
AT  represents  the  change  of  temperature  in  the  concrete.  Substituting  the  linear 
elastic  stress-strain  relationship  in  to  the  linear  momentum  balance  equation,  we  have 


Computation  of  effective  stresses  using  the  principle  of  superposition  described 
above  requires  that  pore  pressures  be  combined  with  stresses  are  induced  by  thermal 
dilation  and  boundary  condition  restraint.  Here,  finite  element  based  thermo- 
mechanical analysis — using  the  ADINA  finite  element  code  (ADINA  2002) — is  used 
for  the  purpose  of  computing  these  thermally  induced  stresses  inside  a circular 
concrete  column  exposed  to  a fire  modeled  using  ASTM  El  19  (ASTM  1995). 

5.5.1  Initial  and  boundary  conditions 

The  initial  conditions  needed  to  solve  for  the  coupled  thermo-mechanical 
equation  of  a concrete  structural  system  are  specified  as  the  absolute  temperature  field 
T at  an  initial  time  / = 0 in  the  domain  Q and  of  boundary  condition  on  the  surface 
f . The  initial  temperature  field  is  then  specified  as 


div  (K-G^S^+lGe'^-lKa^S^+p^S,,  =0 


(5.16) 


5.5  2-D  Axisymmetric  Finite  Element  Stress  Analysis 
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T(x,y,z,Sw,t)[=Q  = T0(x,y,z)  in  Q 


(5.17) 


where  the  initial  temperature  field  T0(x,y,z)  is  a function  of  the  spatial  coordinates 
only.  The  initial  temperature  distribution  T0(x,y,z ) in  Eqn.  (5.17)  is  assigned  with  a 

constant  value  for  the  uniform  temperature,  e.g.,  an  ambient  temperature  of  20  °C. 
The  values  of  temperature  at  the  boundary  Tr  are  specified  to  vary  with  time, 


i.e.. 


T = T(x,y,z,t)  on  Tt 


(5.18) 


The  thermo-mechanical  stress  analyses  conducted  here  are  considered  to  be  one- 
directionally  coupled  because,  as  was  noted  earlier,  temperature  changes  are  assumed 
to  induce  thermal -dilational  deformation  but  mechanical  deformations  are  assumed  to 
produce  no  change  in  system  temperature.  This  assumption  permits  us  to  compute  the 
time-varying  temperature  field  separately,  and  then  to  apply  that  temperature  field  to 
finite  element  thermal  stress  models  to  evaluate  time-varying  thermal  stresses. 

5.5.2  Axisymmetric  finite  element  modeling 

2-D  axisymmetric  simulation  models  are  constructed  for  a 0.3  m diameter 
circular  concrete  column  made  of  a high-strength  concrete  that  is  exposed  to  an 
axially  symmetric  radiant  fire  as  shown  in  Figure  5-1.  TOUGH2  and  the  newly 
proposed  fluid-flow  constitutive  laws  discussed  in  Chapter  4 are  used  to  solve  the 
coupled  heat  and  moisture  flow  equations  and  yield  a complete  time-varying 
description  of  the  temperature  field  throughout  radial  direction  of  the  column.  A one- 
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dimensional  TOUGH2  model  shown  in  Figure  5-1  represents  a typical  radial  slice  of 
the  example  circular  column. 

Computations  are  carried  out  using  the  thermo-mechanical  material  properties 
presented  in  Table  5-1.  Material  parameters  of  high-strength  concrete  under  static 
loading  were  approximated  using  the  assumptions,  1)  the  compressive  strength  of 
plain  concrete  is  approximately  55.0  Mpa  (/c'  = 8000  psi ),  2)  modulus  of  elasticity 
considered  is  for  static  loading  rather  than  dynamic  rates  of  loading,  3)  for  normal- 
weight  concrete  with  a density  of  2200  Kg/m3  ( = 135  lb / ft 3 ),  the  initial  modulus  of 
elasticity  can  be  evaluated  as  E = 57000 . 


Table  5- 

Thermo-mechanical  material  properties 

Material  parameter  and 
description 

Concrete 

Steel 

Units 

E 

(modulus  of  elasticity) 

35000.0 

206850.0 

MPa 

P 

(mass  density) 

2200.0 

7850.0 

Kg/m3 

V 

(Poisson’s  ratio) 

0.19 

0.29 

- 

k 

(heat  conductivity) 

2.9 

46.0 

Watt  l{mC°) 

(heat  capacity) 

921 

419 

Joule  /(kg-  C°) 

a 

(linear  expansion  coeff.) 

11.0E-6 

12.0  E-6 

\IC° 

Temperature  data  from  the  one-dimensional  TOUGH2  models  are  then  mapped 
onto  nodes  located  along  radial  lines  in  the  axisymmetric  finite  element  thermal  stress 
models.  For  each  finite  element  model  nodal  point,  a complete  time-history  of 
temperature  is  extracted  from  the  TOUGH2  simulation  results  and  is  then  applied  to 
the  finite  element  model  as  a time-varying  prescribed  nodal  temperature. 
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Z 


a)  One  dimensional  TOUGH2  model 


z 


150  mm 


z 


c)  Fixed  boundary 
condition  at  top 
of  column 


Z 


d)  Free  boundary 
condition  at  top 
of  column 


b)  Complete  column 


Figure  5-1  A one  dimensional  TOUGH2  model  and  axisymmetric  finite  element 

models  of  concrete  column 
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Time-varying  thermal  stress  analyses  are  then  performed  to  compute  time- 
histories  of  thermal  stresses  throughout  the  two-dimensional  axisymmetric  meshes. 
Thus,  the  axisymmetric  radiant  heat  flux  boundary  condition  is  modeled  as  prescribed 
nodal  temperature  boundary  conditions.  The  mass-  and  energy-balance  equations 
have  been  simultaneously  solved  using  a finite  difference  model  that  simulates 
coupled  one-dimensional  mass  and  heat  transport  phenomena.  Thus,  using  the  time 
history  of  temperature  distributions  obtained  from  the  finite  difference  simulations, 
temporal  nodal  temperatures  of  the  finite  element  model  are  prescribed. 

Due  to  axial  symmetry  of  the  geometry  and  the  thermal  loading  conditions, 
partial  symmetry  finite  element  models  were  used.  In  order  to  investigate  the  effects 
of  boundary  restraint  on  dilatation-induced  thermal  stress  in  the  concrete  column,  two 
different  mechanical  boundary  conditions  are  considered.  The  finite  element  model 
shown  in  Figure  5-lc  shows  axial  constraint  imposed  on  the  both  ends  of  the  column 
and,  along  the  top  boundary,  displacement  in  the  radial  direction  is  constrained  as 
well.  However,  no  constraint  along  the  radial  direction  is  imposed  on  the  model 
shown  in  Figure  5- Id. 

5.5.3  Simulation  results 

Figure  5-2  shows  maximum  principal  stress  variations  sampled  along  Line  A as 
shown  in  Figure  5- lb  and  Figure  5-lc.  Dilatation-induced  thermal  stresses  developed 
in  both  models  increase  as  heating  continues  from  time=300  seconds  to  time=600 
seconds.  More  noticeably,  the  magnitude  of  the  maximum  principal  stresses 
significantly  increases  when  coupling  effects  of  moisture  and  heat  transfer  are  not 
included.  Stress  data  computed  with  consideration  of  moisture  effects  are  generally 
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much  less  severe  than  corresponding  data  computed  without  consideration  of 
moisture  effects. 

A significant  amount  of  thermal  energy  is  consumed  by  latent  heat  during 
evaporation  of  liquid  water  in  concrete  pores.  Phase-transition  and  migration  of 
moisture  during  exposure  to  fire  reduce  the  temperature  increase  inside  partially 
saturated  concrete.  Alternatively,  where  less  severe  thermal  gradients  develop,  lesser 
thermal  gradient  stresses  are  induced.  Both  migration  and  phase-transition  of 
moisture  are  found  to  alter  the  temperature  gradients  developed  particularly  in  near- 
surface zones  of  the  concrete  column. 

Secondly,  comparison  of  the  stress  variation  between  Figure  5-2a  and  Figure  5- 
2b  reveals  that  development  of  the  maximum  principal  stresses  can  change  noticeably 
in  response  to  the  constraints  imposed  on  the  boundaries  of  the  structural  system. 
More  severe  tensile  stresses,  normal  stresses  to  the  axis  of  symmetry  in  the  near 
surface  zone  of  the  column  are  developed  when  the  top  boundary  of  the  column 
model  is  fully  constrained. 

Shown  in  Figure  5-3,  moisture  in  concrete  pores  involved  in  heat  transport  is 
once  again  indicated  by  noticeable  changes  in  the  minimum  principal  stresses.  It  is 
found  that  large  hoop  compressive  stresses  in  near  surface  regions  of  the  column  are 
developed.  The  hoop  stress  component  mainly  contributes  to  significant  development 
of  the  minimum  principal  stresses.  That  is,  maximum  compressive  stress  develops  in 
the  0 -direction  of  the  circular  column.  Furthermore,  the  values  of  the  minimum 
principal  stress  developed  in  the  desaturated  concrete  model  are  much  larger  than 
those  developed  in  the  partially  saturated  concrete  model  (see  Figure  5-3). 
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Constraint  imposed  on  the  top  boundary  of  the  column  is  found  to  have  minimal 
effects  on  the  development  of  compressive  stresses  in  near  surface  regions  of  the 
column.  In  contrast,  this  boundary  condition  has  more  prominent  influence  in  the 
development  of  tensile  stresses.  Temperature  gradients  developed  in  the  radial 
direction  in  the  desaturated  model  are  larger  and  thus,  more  compressive  thermal 
stress  is  induced  in  0 -direction.  It  is  also  noted  that  the  contributions  of  pore  pressure 
to  the  effective  tensile  stresses  are  nearly  negligible. 

Large  tensile  radial-stress  in  the  r-direction  of  the  outer  region  are  developed 
and  found  amplified  when  the  top  boundary  is  restrained  against  vertical 
displacement.  Interestingly,  larger  compressive  hoop-stresses  in  the  9 -direction  of 
the  inner  region  where  a uniform  temperature  rise  occurs  are  developed  regardless  of 
whether  the  boundary  constraint  is  imposed  on  the  top  boundary  of  the  column 
model.  The  outer  shell  of  the  column  thermally  expands  much  more  than  the  cooler 
inner  core  inducing  bending  in  the  column. 

5.6  3-D  Finite  Displacement/Temperature  Mixed  Formulations 

As  discussed  in  Chapter  4,  the  fire  resistance  of  structural  concrete  elements 
such  as  beams  and  columns  can  be  significantly  influenced  by  thermodynamic  states 
in  the  system  during  transient  coupled  moisture  and  heat  transfer  phenomena.  Thus, 
focus  was  placed  primarily  on  quantifying  thermodynamic  state  variables  such  as 
pore-pressure,  temperature,  and  saturation  that  characterize  the  hydro-thermal 
behavior  of  concrete  under  a thermal  loading  condition.  In  addition,  the  effects  of 
steel  reinforcement  of  the  structural  concrete  elements  on  hydro-thermal  behavior  of 
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a)  Maximum  principal  stresses  (fixed  boundary) 


b)  Maximum  principal  stresses  (free  boundary) 

Figure  5-2  Comparison  of  maximum  principal  stresses  at  t = 600  sec. 


Minimum  principal  stress  (MPa)  Minimum  principal  stress  (MPa) 
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a)  Minimum  principal  stresses  (fixed  boundary) 


a)  Minimum  principal  stresses  (free  boundary) 

Figure  5-3  Comparison  of  minimum  principal  stress  variations  with  respect  to  initial 

and  boundary  conditions 
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concrete  under  a thermal  loading  condition.  In  addition,  the  effects  of  steel 
reinforcement  of  the  structural  concrete  elements  on  hydro-thermal  behavior  of 
concrete  were  investigated  and  quantitatively  identified  for  the  influence  of 
heterogeneity  of  the  system  in  the  development  of  thermodynamic  states. 

In  subsequent  sections,  a three-dimensional  finite  element  with  temperature 
degree  of  freedom  will  be  developed.  Transient  temperature  profiles  that  were 
obtained  from  three-dimensional  finite  difference  model  of  simulating  hydro-thermal 
behavior  of  reinforced  high-strength  concrete  column  exposed  to  a hydrocarbon  pool 
fire  in  Chapter  4 is  then  mapped  into  the  finite  element  domain  to  perform  thermal 
stress  analysis. 

Recall  the  volume-averaged  linear  momentum  balance  equation,  Eqn.  (5.16). 
For  simplicity  in  the  following  derivation,  the  gravitational  effects  term  is  neglected. 
To  solve  the  equation,  employ  a weighted  residual  method  by  introducing  residual 
Rn  due  to  the  error  in  the  approximation 

<519 

Introducing  a weak  formulation  by  using  a weighting  function  wt  = w(  (x) , we  can 
express  the  weighted  residual  as 


(5.20) 


Thus,  we  have 


[wdUjdn  =o 


(5.21) 
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Using  integration  by  parts,  the  weighted  residual  equation  becomes 


{ (W'dj) jdQ—  [ 


WijdjdD.  = 0 


(5-22) 


By  applying  the  divergence  theorem  to  the  above  equation,  it  is  rewritten 


Wij^cin  = o 


(5.23) 


or 


£kedQ- | wi  j(-3Ka)SiJ(Af)dQ.  = 0 


(5-24) 


where  £ and  T are  trial  functions  that  satisfy  the  essential  boundary  conditions  on 
the  boundary  of  the  system,  the  stress  vector  is  denoted  by  t,  = djn] , nj  is  a surface 

normal  to  the  boundary  surface  T , K represents  the  bulk  modulus,  and  a thermal 
expansion  coefficient  is  denoted  by  a . This  is  a weak  form  of  the  linear  momentum 
balance  equation  for  a coupled  thermo-elastic  problem. 

Discretization  of  the  above  equation  is  performed  by  using  the  approximated 
solution  uh  as  the  same  weighting  function  w 

h (5.25) 

w = u ' 

that  must  satisfy  the  essential  boundary  conditions.  The  approximated  solution  uh  can 
then  be  represented  in  terms  of  nodal  displacements  d:  and  nodal  temperatures  Tt  as 
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we  perform  numerical  integration  using  Gauss  quadrature  over  the  element  volume 
in  the  natural  coordinate  system,  i.e.,  using  local  axes  r,  s , and  t,  which  results  in 


[/r,]{d}+[Mr]{rJ  = {F„| 


(5.32) 


where  [£,,]  is  the  mechanical  stiffness  matrix,  [d]  denotes  the  global  mechanical 


displacement  vector,  [ Mr ] represents  the  thermo-mechanical  coupling  matrix,  {ro} 

denotes  the  global  nodal  temperature  vector,  and  {Fexlemal}  represents  the  surface 

force  vector.  More  explicitly,  numerical  integration  is  used  to  evaluate  the  multi- 
dimensional integrals  in  Eqn.  (5.32)  as  follows 

NT 


[*.]  = ! f[i*nc][BpQ 
= Z ['SX^HBUcWdrdsd, 


e=l 


NT 


=z 


e=\ 


/=l,int  7=1, int  k=\, int 

Z Z £(Sr(r„»JA)'C'S(r„J),ll)-det)-W,.Wi.W1) 
a a n. 


(5.33) 


= Z f,'£lf,ldetJ[fifl fiWdrdsdt 


e=l 
NT 

-z 

e=\ 


/=l,int  7=1, int  /c =1  ,int 

Z Z Z ■/9-Af(ri,sJ.l,)detJ-W,.Wy.W<) 

q,  nt  a. 


(5.34) 


where  NT  represents  the  total  number  of  elements,  “ det  J ” is  the  determinant  of  the 
Jacobian  matrix  transforming  global  to  natural  coordinates  (Bathe  1 996),  the  number 
of  integration  sampling  points  in  each  element  is  denoted  by  “int”,  and  the  W’s 
represent  the  weights  corresponding  to  the  Gauss  Quadrature  sampling  points.  The 
surface  force  vector  {Fexlemal}  consists  of  the  following  three  parts 
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and 


ra=£  [m  wrfr, 

e=l 

NT 
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(5.37) 
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i=int  7=int 
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where  AT  represents  the  total  number  of  elements,  det  Jy  is  the  determinant  of  the  2- 

D Jacobian  matrix  transforming  global  to  natural  coordinates  on  element  boundaries 
(/  = x,y,z)  (Bathe  1996). 


5.7  Discretization  of  the  Transient  Heat  Conduction  Equation 

Having  established  the  Galerkin  form  of  the  weighted  residual  method  when 
applied  to  the  linear  momentum  balance  equation,  we  now  employ  the  same  method 
of  weighted  residuals  to  form  a weak  formulation  of  the  volume-averaged  transient 
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heat  conduction  equation  given  in  Eqn.  (5.3).  However,  in  order  to  solve  the  partial 
differential  equation  of  heat  conduction,  we  need  the  specification  of  the  initial 
condition  at  time  t =t0  in  initial  domain  Q and  the  specification  of  the  boundary 

condition  on  the  surface  T for  the  proposed  problem.  The  initial  temperature  field 
considered  in  this  study  is  specified  as 


T(x,y,z,0)  = T0(x,y,z)  = 0 in  Q 


(5.38) 


Also,  the  values  of  temperature  at  the  boundary  Tr  are  specified,  which  may  be 
constant  or  may  be  allowed  to  vary  with  time 


T = T(x,y,t)  in  Tt 

and  which  are  known  as  a Dirichlet  or  essential  boundary  condition.  In  the  other 
condition,  the  values  of  the  heat  flow  in  the  direction  n normal  to  the  boundary  T q are 

prescribed  as  q B(x,y,z,t) 


, o dT 

q8=~kS,J!h  on 


(5.40) 


which  is  called  a Neumann  or  natural  boundary  condition. 

The  problem  of  heat  transfer  in  a concrete  structure  exposed  to  fire — where  the 
temperature  field  at  various  points  in  the  domain  varies  with  time — is  considered  to 
be  a transient  problem.  For  such  problems  a temporal  discretization  for  the  time 
derivative  of  the  transient  heat  conduction  is  required  in  addition  to  the  spatial 
discretization.  In  the  following  section,  temporal  discretization  methods  are  presented 
for  the  finite  element  analysis  of  transient  heat  conduction  problems. 
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The  residual  for  the  energy  balance  equation  is  defined  as 


da, , 

Rn=(kT,\,+Q  + T ,J 


dT 

dT  eu  PCp  dt 


(5.41) 


Prior  to  a further  development  of  temporal  discretization,  it  is  once  again  noted  that 
the  concrete  structure  considered  in  this  study  can  be  reasonably  analyzed  without 


having  to  consider  the  coupling  factor 


da  a 

rj,  IJ 


dT  ,J 


in  the  residual  for  the  energy 


balance  equation.  Thus,  we  have 


Ra  =(kfl)l+Q-pC  — 

dt 


(5.42) 


The  subsequent  steps  of  the  derivation  to  obtain  the  Galerkin  form  of  weighted 
residual  expression  remain  same  as  outlined  in  the  previous  section.  The  only 
exception  in  the  weak  form  of  the  transient  heat  conduction  is  the  temporal  derivative 


term.  As  defined  in  Eqn.  (5.26),  a set  of  linear  independent  functions  [iV]  is  used  as 

finite  element  basis  functions  in  the  strategy  of  dividing  the  solution  domain  into 
finite  elements  and  using  the  element  shape  functions  to  define  the  solution  within  the 
element  that  is  based  on  the  nodal  values  of  temperature.  In  this  case,  however,  the 
shape  functions  are  also  applied  to  the  temporal  derivative,  i.e.. 


dTe  . dT‘  * 

~—(x,y,z)=  2,  -T-N,(x,y,z) 
dt  a dt 


(5.43) 


137 


where  NP  represents  the  number  of  nodes  in  each  element  and  fe  represents 
temperature  residual  in  each  element  occupying  an  element  domain  Qt, . This  results 
in  the  following  Galerkin  finite  element  form  of  the  discretized  equations 


[{B]T[k][B]dn{T0}+l[N}T pCr[N]dn{Ta}=l'  [N^{qB}dV  <5'44) 

where  [&]  represent  the  thermal  conductivity  coefficient  matrix,  [5],  the  gradients 

of  the  basis  functions  as  defined  in  Eqn.  (5.28),  and  the  heat  flux  on  the  surface 
boundary  T is  defined  using  the  surface  normal  n as 


,,  dT 

q,=-ks“!fo  °"  r» 


(5.45) 


Equation  (5.44)  represents  the  spatially  discretized  finite  element  form  of  the  transient 
heat  conduction  equation.  Thus,  in  a more  compact  form,  we  rewrite  Eqn.  (5.44) 


(5.46) 


where  [ CH  ] represents  the  global  heat  capacity  matrix,  [/C7  ] is  the  global  heat 
conductivity  matrix,  and  {//„,}  represents  the  global  thermal-loading  vector  after  the 
multi-dimensional  integrals  are  numerically  evaluated  as 

NT 

e=\  e 

= t f'f'f (5.47) 


e=\ 
NT 

-z 

e=\ 


/=l,int  y=l,int  Ar=l,int 

a n.  n,  ’ 


138 


and 
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The  global  thermal-load  matrix  { } consists  of  the  following  three  parts 
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The  mapping  functions  [/V]  used  in  this  study  for  the  nodal  temperature  degrees  of 
freedom  are  first  order  polynomial  functions.  Two-point  Gauss  quadrature  is  then 
used  to  integrate  the  element  heat  capacity  matrix  [ CH  ] and  the  element  heat 
conductivity  matrix  [ATr] . 


5.8  27/8  Finite  Element 

As  was  mentioned  earlier,  lower-order  temperature  interpolations,  i.e.,  first 
order  polynomial  functions,  have  been  used  as  the  finite  basis  functions  for  nodal 
temperature  degrees  of  freedom,  but  higher-order  displacement  interpolations,  i.e., 
second  order  polynomial  functions,  have  been  used  for  nodal  displacement  degrees  of 
freedom.  Discussion  of  appropriate  temperature  and  displacement  interpolations  that 
should  be  used  in  order  to  develop  an  effective  element  follows. 

Assume,  for  instance,  that  we  use  27  nodes  to  interpolate  the  displacements  in  a 
three-dimensional  element  and  we  assume  linear  variations  of  the  temperature  field 
by  using  only  8-node  temperature  interpolation  within  the  same  element.  Such  an 
element  is  referred  to  as  a “27/8  element”.  In  this  formulation,  the  element 
temperature  is  defined  by  nodal  temperature  variables  that  pertain  to  adjacent 
elements  in  the  assemblage.  Thus,  continuity  of  temperature  between  elements  is 
always  contained  in  the  solution.  Then,  the  temperature  prediction  could  be  of  higher 
order  and  could  therefore  hypothetically  be  more  accurate. 

However,  studies  have  shown  (e.g.,  Bathe  1996)  that  when  a mixed  finite 
discretization  such  as  the  u/p-c  formulation  is  used,  then  a field  variable  (e.g.,  pore 
pressure,  temperature)  coupled  with  mechanical  displacement  should  not  be 
interpolated  at  too  high  a degree  because  locking  could  be  introduced  into  the  finite 
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element.  Thus,  higher-order  interpolation  for  the  field  variable  associated  with 
coupling  of  the  mechanical  displacement  field  would  have  a negative  effect  on  finite 
element  prediction  of  the  mechanical  displacement  field. 

Secondly,  the  thermo-mechanical  coupling  matrix  [ M ] that  is  numerically 
integrated  for  the  product  of  temperature  and  mechanical  displacement  interpolation 
functions  would  be  overly  stiffened  when  Gauss  points  of  a higher-order  rule  are  used 
to  increase  the  accuracy  of  a finite  element  analysis.  The  number  of  Gauss  points  has 
a lower  limit  because  element  volume  has  to  be  integrated  exactly  in  the  limit  of  mesh 
refinement.  Thus,  when  exact  integration  is  used,  then  small  change  in  the 
temperature  field  combined  with  a higher-order  integrated  coupling  matrix  can 
numerically  induce  large  internal  force.  As  a result,  despite  the  use  of  an  exact 
integration  scheme  to  form  the  thermo-mechanical  coupling  matrix  [M] , the  finite 
element  solution  of  the  coupled  thermo-mechanical  stress  field  may  not  be  reliable. 
This  can  occur  if  over-estimation  of  the  mechanical  displacement  field  results  in  a 
large  internal  force  computation  (Cook  et  al.  1989). 

Therefore,  we  want  to  use  the  highest  degree  of  temperature  interpolation  that 
neither  introduces  locking  into  the  element  nor  causes  poor  quality  of  the  temperature 
prediction.  The  27/8  element  shown  as  Figure  5-4  is  most  suitable  for  this  purpose 
because  (1)  its  mixed  finite  element  discretization  yields  a reliable  solution,  i.e.,  it  is 
stable  and  convergent  (Bathe  1 996)  and  (2)  the  rate  of  convergence  of  the 
temperature  field  (and  thus  thermally  induced  stresses)  as  the  mesh  is  refined  is  of 
order  0(h 2)  (Huang  and  Usmani  1994). 
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• Node  with  displacement  variables 
® Node  with  displacement  and  temperature  variables 

Figure  5-4  A schematic  sketch  of  effective  27/8  element 

5.9  Solution  Approach 

In  this  section,  the  spatially  discretized  finite  element  formulation  of  the 
transient  heat  conduction  equation — a first  order  system  of  ordinary  differential 
equations  with  respect  to  time — will  be  also  discretized  in  the  time  domain.  The  finite 
difference  discretization  scheme  in  which  the  time  domain  is  divided  into  small  time 
increments  is  used  to  solve  the  equation  step  by  step  in  the  given  space  domain.  The 
objective  is  to  develop  an  algorithm  where  both  the  space  and  time  domains  are 
discretized  in  a piecewise  manner.  Thus,  the  unknown  values  dn+l  and  Tn+i  in  Eqns. 


(5.32)  and  (5.44)  at  time  tn+i  will  be  determined  where  the  values  Tn  are  known  at 
time  tn  and  so  is  Hexl  in  the  time  interval  A / . 


142 


First,  we  combine  Eqn.  (5.32)  with  Eqn.  (5.44)  in  matrix  form 
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where  three  displacement  degrees  of  freedom  are  allowed  at  all  27  nodes  of  the  27/8 
finite  element  (shown  in  Figure  5-4),  but  only  8 comer  nodes  have  temperature 
degree  of  freedom.  In  alternate  form,  this  may  be  expressed  as 
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In  subsequent  discussions,  it  is  convenient  to  condense  Eqn.  (5.53)  into  the  following 
form 
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Therefore,  it  becomes  a first  order  ordinary  differential  equation  (ODE)  and  we  now 
determine  the  solution  of  an  initial  value  problem  Given  d0  = d{ 0) , find  d = d(t) , 

such  that  A d+  g(d ) = h(t) . Eqn.  (5.54)  is  a system  of  first  order  ODEs  that  can  be 
numerically  solved  by  discretizing  the  time  domain  into  finite  time  intervals, 

A tk  = tk+]  -tk , where  the  temporal  solutions  are  approximated  as  dk+]  ~ d(tk+x ) and 

Defining 

d'  = Ad  , d"  - Ad  , 

and  (5-57) 

/*  =AAd,t)  = h(t)-g(d) 

we  transform  Eqn.  (5.54)  into 

d * = f\d,t ) (5-58) 

Introducing  the  general  linear  k-step  method  (Borja  1991)  that  approximates  the 
solution  of  ODEs  of  the  first  order 

k 

Z (amdl i-„  + At ) = o (5.59) 

m= 0 

where  am  and  pm  are  unknown  coefficients.  Therefore,  substituting  Eqn.  (5.57) 
back  into  the  above  equation  gives 

k k 

^Z  K<+.-m)  + Z At  P„XK+\-m  ~Sn+l-m)  = 0 (5.60) 

01=0  m=0 


which  is  equivalent  to  the  following  matrix  form 
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Also,  in  component  form,  Eqn.  (5.60)  is  equivalent  to 

[*;.]{‘7L,-.+[wrl{?;L,.={F„} 


n+\-m 


(5.62) 


and  letting  a0  = -1  and  /?0  = /?,  = /?2  = = J3k  and  J30*  0 , we  have 

[C„]f  {?;}„„  - £ a „ 1+ A/A  = a/ a {//„, } 

V m= 1 2 


n+1 


(5.63) 


Because  the  equilibrium  components  of  (hn+l_m  - gn+1_m)  for  m=  1 , 2,. . .,k  vanish 
from  the  previous  time  step  tn , the  linear  momentum  and  the  energy  balance  must  be 
ensured  at  {^}n+l  such  that 


/(4m)  = 


I n+l 


{A»L 

a/a  -a/a  K]{r,}„,  -a 


= 0 (5.64) 


where  V(\TA,„)  = ~\  \TA,.rY.  ]■  Etln'  (5-64) isthen  solved  by 

A>  V 


m= 1 


using  the  Newton- Raphson  method,  linearizing  Eqn.  (5.64)  to  solve  for  dJn++\  of  the 
y+ 1-th  iteration  yields 


A//.. 


m:\ ) - /«, ) + (a;;-a„) + o<*!  > 


ad. 


»+l 


(5.65) 


When  it  converges,  that  is,  / (//,')  « 0 , the  formula  becomes 


A/+i  = 


g/(^+1)V' 

V 7 


/(O 


(5.66) 


Or, 
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MK.H'i,} 


(5.67) 


where  (A'**]  represents  a global  stiffness  matrix  that  contains 


'WjJ' 

V ^dn+\  j 


Kl 


[W] 


[0]  A tPa[KT]+[CH] 


(5.68) 


and  the  residual  force  vector  {r}  is  defined  as 


L 


{FJ„AK.MnAMr\{TX 


(5.69) 


Thus,  if  | ,f  at  the  current y-th  iteration  does  not  satisfy  the  equilibrium,  e.g., 

Jr,/, | j * 0 or  Jr^|  j > tolerance , then  it  will  be  updated  as  \dj~l  j - Jef,/,,  j + { .V/,.,  j in 
subsequent  iteration  until  convergence  is  reached. 


5.10  3-D  Finite  Element  Thermoelastic  Stress  Analysis 

In  this  section,  transient  thermo-mechanical  stress  analysis  is  performed  using  a 
three-dimensional  (3-D)  finite  element,  i.e.,  the  27/8  element  shown  in  Figure  5-4, 
which  has  been  developed  and  coded  in  the  PLASFEM.  In  the  following,  a reinforced 
concrete  column  exposed  to  fire  will  be  modeled  using  this  3-D  finite  element. 

5.10.1  Model  description 

A square  cross-section  (380x380  mm)  reinforced  concrete  column  exposed  to  a 
radiant  fire  condition  shown  in  Figure  5-5  (see  also  Chapter  4)  is  modeled.  The  size 
of  longitudinal  bars  is  a no.  6 bar  and  ties  of  no.  3 bars  are  provided  and  served  as 
shear  reinforcement  for  the  reinforced  concrete  column.  The  vertical  spacing  of  ties 
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Figure  5-5  A reinforced  concrete  column  exposed  to  radiant  heat 


used  in  the  model  is  40  mm.  The  ends  of  the  ties  are  lap  sliced  with  comer  hooks  that 
are  idealized  as  perfectly  anchored  by  a sharp  90°  bend  around  a bar.  Confining 
effects  in  the  concrete  core  and  interface  between  the  longitudinal  bars  and  the  ties 
are  not  considered  in  the  finite  element  analysis.  More  importantly,  the  effects  of  steel 
reinforcement  on  hydro-thermal  behavior  of  concrete  are  considered  in  order  to 
quantitatively  identify  the  influence  of  heterogeneity  of  reinforced  concrete  structural 
elements  on  the  thermodynamic  state  of  the  system  and  thus,  the  influence  on  the 
development  of  thermally  induced  stress  field. 

Secondly,  because  of  the  symmetry  of  the  geometry  and  the  thermal  loading 
conditions,  a partial  symmetry  model  of  the  actual  column  segment  is  developed 
(illustrated  in  Figure  5-6)  as  a single  quadrant  mesh  symmetric  in  both  the  xz-  and  the 
xy-plane.  Figures  5-6a  and  5-6c  illustrate  the  finite  element  mesh  layout  in  a plane 
view  at  two  different  y-coordinates,  i.e.,  y=2.5  mm  and  y=20.0  mm,  respectively.  In 


subsequent  discussion,  the  plane  view  shown  in  Figure  5-6a  is  referred  as  Plane  A 
while  the  one  shown  in  Figure  5-6b  is  referred  as  Plane  B. 
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(c)  (d) 

Figure  5-6  3-D  finite  element  models  of  the  concrete  column  segment 


As  discussed  earlier  in  this  chapter,  the  degree  of  constraint  provided  by  the 
boundaries  of  the  column  segment  will  affect  the  thermal-dilation  related  stresses  that 
are  computed.  Once  again,  two  different  sets  of  boundary  conditions  that  bracket 
realistic  field  conditions  (e.g.,  material  constraint)  are  considered.  In  the  first  case 
(Figure  5-6c),  vertical  expansion  is  completely  restrained,  while  in  the  second  (Figure 
5-6d),  the  top  of  the  column  segment  is  free  to  expand.  At  a plane  of  geometric 
symmetry,  mechanical  displacement  boundary  conditions  for  the  symmetric  thermal 
loading  are  modeled  such  that  no  translational  motion  is  allowed  perpendicular  to  a 
plane  of  geometric  symmetry  (Figure  5-6). 
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5.10.2  Modeling  of  the  transient  thermal  loading  condition 

Using  TOUGH2  to  solve  the  coupled  heat  and  moisture  flow  equations  yields  a 
complete  time-varying  description  of  the  temperature  field  throughout  three- 
dimensional  directions  of  the  column.  Temperature  data  from  the  three-dimensional 
TOUGH2  model  (refer  to  Chapter  4)  is  then  mapped  onto  nodes  located  in  the  finite 
element  thermal  stress  models.  For  each  finite  element  model  nodal  point,  a complete 
time-history  of  temperature  is  extracted  from  the  TOUGH2  simulation  results  and  is 
then  applied  to  the  finite  element  model  as  a time-varying  prescribed  nodal 
temperature  {r/re4c"w} . Time-varying  thermal  stress  analyses  are  then  performed  for 

each  boundary  condition  model  to  compute  time-histories  of  thermal  stresses 
throughout  the  three-dimensional  meshes.  At  each  time  step,  the  following  system 
equation  is  solved 

={frL  (5.70) 

where 

{%„  = (5.71) 

Therefore,  we  compute  thermal  stresses  at  time  tn+l  through  the  relationships 

{°}n^=[C]{s}n-{P}  ATn+1  (5.72) 

where  {e}n  = [B}{d]n  and  ATn+l , temperature  change  at  a point  of  interest  in  the 
natural  coordinate  system  at  time  tn+l , is  computed 


(5.73) 
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5.10.3  Simulation  results 

As  demonstrated  earlier  in  the  case  of  2-D  axisymmetric  finite  element  stress 
analysis,  computation  of  effective  stresses  using  the  principle  of  superposition 
described  by  Eqn.  (5.12)  requires  that  pore  pressures  be  combined  with  stresses  that 
are  induced  by  thermal  dilatation  and  boundary  condition  restraint.  First,  the  thermal- 
dilatation  total  stress  computed  by  finite  element  analysis  is  taken  as  cr  since  it 
represents  the  macroscopically  averaged  effect  of  thermal  dilation  and  boundary 
constraint.  Then,  Eqn.  (5.12)  is  used  to  combine  (or  superimpose)  thermal  stresses 
with  pore  pressures  separately  computed  using  the  3-D  finite  difference  model 
described  in  Chapter  4.  Thus,  in  the  subsequent  discussion,  both  maximum  and 
minimum  principal  stresses  represent  principal  values  of  the  effective  stress  tensor 
a" . 

Figure  5-7  and  Figure  5-8  show  maximum  principal  stresses  cr,  developed  in 
the  case  of  a fully  constrained  boundary  condition  imposed  on  y-direction 
displacements  shown  in  Figures  5-6c.  This  case  will  be  referred  as  the  fixed  boundary 
condition  case.  Maximum  principal  stresses  cr,  computed  in  both  Plane  A (Figure  5- 
6a)  and  Plane  B (Figure  5-6b)  vary  in  a range  from  40  MPa  to  80  MPa.  at  time  t=200 
sec.  after  all  the  four  sides  of  the  column  initially  exposed  to  a fire.  Owing  to  large 
effective  tensile  stresses  acting  normal  to  the  xz-plane,  i.e.,  cr”y , the  severe  maximum 

principal  stresses  acting  on  the  principal  xz-plane  are  developed  in  the  concrete  cover 
zone  of  the  column  segment.  One  particularly  noteworthy  feature  of  the  thermally- 
induced  stress  field  is  severe  development  of  effective  tensile  stress  cr",  in  the  near- 
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surface  regions  along  the  surface  of  the  column  segment  and  in  the  near-comer 
regions. 

Figures  5-9  and  5-10  show  maximum  principal  stresses  cr,  developed  in  the 
case  of  no  constraint  imposed  on  the  top  boundary  of  the  column  segment  as  shown 
in  Figures  5-6d.  In  contrast  to  the  fixed  boundary  condition  case,  maximum  tensile 
stress  components  of  the  effective  stress  tensor  computed  in  the  near-surface  regions 
are  cr”t  and  cr”  rather  than  a”y . Thus,  large  effective  tensile  stresses  are  developed 

perpendicular  to  the  yz-plane,  i.e.,  cr"a,  and  mainly  contribute  to  the  maximum 
principal  stress  in  the  surface  region  along  the  x-coordinate.  However,  acting 
perpendicular  to  the  xy-plane,  cr”  is  found  to  be  the  dominant  tensile  stress 
component  along  the  z-coordinate.  This  reveals  that  characteristics  of  the  stress  field 
developed  in  the  free  boundary  condition  case,  the  development  of  dominant  tensile 
stress  component  is  direction-dependent,  but  results  in  the  development  of  large 
maximum  principal  stresses  in  the  surface  regions.  It  is  noted  that  the  peak  value  of 
maximum  principal  stress  in  the  comer  region  is  only  one-quarter  the  magnitude  of 
what  is  computed  in  the  fixed  boundary  condition  case.  Thus,  in  the  fixed  boundary 
condition  case,  more  severe  tensile  effective  stresses  cr”  is  thermally  induced  in  the 

y-direction  perpendicular  to  the  symmetric  xz-plane  near  the  comer  region  of  the 
column  segment. 

Comparing  these  maximum  principal  stresses  within  the  first  1 0 mm  of  depth  in 
both  boundary  condition  cases,  significant  changes  in  magnitude  and  trend  are  found 
in  the  stress  field.  It  appears  that  less  severe  thermal  stresses  develop  in  the  free 
boundary  condition  case.  No  significant  tensile  stress  was  developed  in  regions  near 
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the  locations  of  the  longitudinal  bars  in  the  free  boundary  condition  case  due  to  free 
expansion  (e.g.,  Figure  5-8  vs.  Figure  5-10).  One  noteworthy  feature  of  thermally- 
induced  stress  fields  computed  for  both  the  boundary  condition  cases  is  that  the 
maximum  principal  values  of  the  stresses  would  exceed  the  ultimate  tensile  strength 
of  a typical  high-strength  concrete — not  accounted  for  in  the  linear  elastic  analyses 
conducted  here — and  failure  would  thus  occur  before  these  stress  levels  could  be 
reached. 

The  influence  of  non-uniform  temperature  field  on  the  development  of 
compressive  stress  is  discussed.  As  already  seen  in  the  2-D  axisymmetric  study,  the 
general  trends  of  stress  contour  of  the  minimum  principal  stress  <x3  at  time  t=200  sec. 

are  found  to  be  quite  similar  in  both  the  fixed  (Figures  5-1 1 and  5-12)  and  the  free 
(Figures  5-13  and  5-14)  boundary  condition  cases.  The  magnitudes  of  compressive 
normal  stresses  cr*t  and  <j".  computed  in  the  concrete  region  in  the  vicinity  of  shear 
reinforcement  (steel  shear  ties)  are  not  significantly  different  in  the  two  cases.  With 
steep  temperature  gradients,  the  outer-surface  regions  expand  and  such  an  expansion 
cannot  proceed  freely  in  the  relatively  cooler  regions,  the  inward  thermal  expansion 
of  the  outer  surface  region  is  suppressed  in  the  xz-plane  by  the  cooler  inner  concrete. 
This  results  in  significant  development  of  compressive  stresses  a[x  and  a”„  in  the 

concrete  in  the  vicinity  of  shear  reinforcement,  which  are  acting  normal  to  the  surface 
exposed  to  fire. 

However,  cr^  that  is  developed  in  the  same  region  of  the  fixed  boundary 

condition  case  is  the  dominant  tensile  stress  component  of  the  effective  stress  tensor, 
which  contributes  to  the  development  of  large  maximum  principal  stresses  cr,  shown 
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(a)  Maximum  principal  stresses  in  Plane  A of  the  fixed  boundary  case 
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(b)  Stress  contour  in  the  Plane  A 

Figure  5-7  Maximum  principal  stresses  at  t=200  sec.  (fixed  boundary) 
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(a)  Maximum  principal  stresses  in  Plane  B of  the  fixed  boundary  case 
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(b)  Stress  contour  in  the  Plane  B 

Figure  5-8  Maximum  principal  stresses  at  t=200  sec.  (fixed  boundary) 
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Maximum  principal  stress 


(a)  Maximum  principal  stresses  in  Plane  A of  the  free  boundary  case 
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(b)  Stress  contour  in  the  Plane  A 


Figure  5-9  Maximum  principal  stresses  at  t=200  sec.  (free  boundary) 
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Maximum  principal  stress 


(a)  Maximum  principal  stresses  in  Plane  B of  the  free  boundary  case 
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(b)  Stress  contour  in  the  Plane  B 

Figure  5-10  Maximum  principal  stresses  at  t=200  sec.  (free  boundary) 
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in  Figure  5-8.  Thus,  a state  of  stress  where  cr,  and  cr3  have  opposite  signs  is 
developed  particularly  in  the  concrete  regions  adjacent  to  shear  reinforcement. 
Subsequently,  the  shear  stress  component  cr", , having  a magnitude  as  significant  as 

tensile  normal  stress  component  a”  is  observed  in  the  effective  stress  tensor 
developed  in  that  region. 

For  the  free  boundary  condition  case,  tensile  cr"  is  almost  negligible  in  the 
near-surface  region  where  the  dominant  tensile  stresses  are  rather  cr"t  and  cr"t . Small 
compressive  cr^  is  developed  in  the  comer  region  where  cr"t  and  cr"t  also  become 
compressive.  Since  three  normal  stress  components  are  compressive  in  the  comer 
region,  a state  of  stress  is  developed  where  maximum  principal  stress  dj  and 
minimum  principal  stress  <x3  have  same  negative  signs.  Thus,  there  is  a small  region 

of  the  inner  comer  where  the  pure  compression  occurs  in  the  free  boundary  condition 
case. 

5.10.4  Summary 

The  boundary  constraint  imposed  on  the  y-direction  displacement  of  the  column 
segment  is  found  to  be  very  influential  in  development  of  a state  of  thermal-dilatation 
stresses  in  the  3-D  model.  The  minimal  differences  in  magnitude  and  the  similarity  in 
trends  between  the  results  of  minimum  principal  stresses  are  found  while  significant 
differences  in  development  of  tensile  stresses  are  noticed  by  observation  in  the 
distributions  of  the  maximum  principal  stresses.  Examination  of  the  effective  stress 
distributions  at  various  spatial-temporal  locations  in  Plane  B in  both  the  fixed  and 
free  boundary  conditions  case  reveals  the  effects  of  the  boundary  condition  on  the 
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(a)  Minimum  principal  stresses  in  Plane  A of  the  fixed  boundary  case 
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(b)  Stress  contour  in  the  Plane  A 

Figure  5-1 1 Minimum  principal  stresses  at  t=200  sec.  (fixed  boundary) 
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(a)  Minimum  principal  stresses  in  Plane  B of  the  fixed  boundary  case 
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Figure  5-12  Minimum  principal  stresses  at  t=200  sec.  (fixed  boundary) 
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(a)  Minimum  principal  stresses  in  Plane  A of  the  free  boundary  case 
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Figure  5-13  Minimum  principal  stresses  at  t=200  sec.  (free  boundary) 
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(a)  Minimum  principal  stresses  in  Plane  B of  the  free  boundary  case 
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(b)  Stress  contour  in  the  Plane  B 

Figure  5-14  Minimum  principal  stresses  at  t=200  sec.  (free  boundary) 
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development  of  the  effective  stresses,  particularly  two  effective  stress  components, 
i.e.,  a"y  and  cr"  . 

First,  the  stress  component  <r”y  in  the  effective  stress  tensor  is  noticeably 

affected  by  the  boundary  constraint.  Steep  temperature  gradient  induced  by  a high 
rate  of  heating  combined  with  the  boundary  restraint  and  internal  pore-pressure 
buildup,  the  effective  stress  component  a"iy , which  is  acting  parallel  to  the  surfaces  of 

fire  exposure,  is  found  to  exceed  a typical  tensile  strength  of  high-strength  concrete  at 
early  stages  during  fire  exposure.  In  addition,  cryy  is  the  dominant  tensile  stress 

component  in  the  development  of  maximum  principal  stress  than  any  other  stress 
components  in  the  fixed  boundary  condition  case. 

Secondly,  the  same  normal  stress  component  & ^ in  the  free  boundary  condition 

case  is  rather  insignificant  compared  to  the  other  normal  stress  components  in  the 
development  of  thermal  effective  stresses.  In  planes  parallel  to  the  surfaces  exposed 
to  fire,  cr”x  acting  on  the  yz-plane  and  cr"  acting  on  the  xy-plane  are  the  dominant 

tensile  stress  components  that  mainly  contribute  to  the  maximum  principal  stresses  in 
the  near-surface  regions.  However,  we  found  these  normal  stresses  less  significant  in 
inner  layers  of  the  concrete  cover.  The  shear  stress  cr  ",  increases  as  those  normal 
stresses  decrease  in  the  vicinity  of  the  shear  reinforcement. 

Therefore,  considering  the  degree  of  restraint  provided  by  the  column 
boundaries,  an  actual  state  of  stress  field  in  realistic  field  conditions  is  thought  to  be 
bracketed  by  the  thermal-dilation  related  stresses  that  are  computed  for  two  different 
sets  of  boundary  conditions.  Hypothetically,  if  brittle  material  failure  of  high-strength 
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concrete  associated  with  thermally  induced  stress  would  take  places  in  the  concrete 
column,  then  it  would  be  tensile  fracture  due  to  severe  development  of  effective 
tensile  stresses,  , normal  to  the  xz-plane  where  fracture  planes  could  be  formed 


most  likely  first  in  the  near-surface  comer  regions  of  the  column. 


CHAPTER  6 

CONTRIBUTION  OF  PORE  PRESSURE  AND  THERMAL  GRADIENT  STRESS  TO 

SPALLING 

Inconsistent  occurrence  of  explosive  spalling  has  been  observed  in  the  past  tests  of 
high-strength  concrete  (Phan  and  Carino  2002).  In  regard  to  the  thermal  behavior  of  high- 
strength  concrete,  the  effect  of  low  permeability  on  transport  phenomena  leading  to  pore 
pressure  buildup  was  suspected  to  be  main  cause  for  higher  susceptibility  of  high-strength 
concrete  to  explosive  spalling.  This  chapter  addresses  the  contribution  of  high  pore- 
pressure  buildup  and  thermal  gradient  stress  to  thermal  spalling  of  high  strength  concrete, 
based  on  the  results  obtained  in  this  study. 

6.1  Factors  Influencing  Spalling  at  Early  Stages  of  Heating 

Based  on  the  linear  thermo-elastic  stress  results  presented  in  Chapter  5,  the  pore 
pressure  associated  with  moisture  transport  phenomena  at  early  stages  of  heating  (during 
the  first  five  minute  fire  exposure)  is  found  to  have  nearly  negligible  effects  on  the 
development  of  the  effective  stress  tensor  compared  to  the  magnitude  of  thermal- 
dilatation  stress  components.  A most  significant  factor  of  spalling  at  an  early  stage  of 
heating  is  thermal  gradient  stresses  that  can  cause  fracture  failure. 

The  effects  of  moisture  on  the  temporal  temperature  distributions  result  in  very 
significant  changes  of  thermally  induced  stress  field.  As  investigated  using  an 
axisymmetric  finite  element  analysis  in  the  previous  chapter,  the  magnitudes  of  thermally 
induced  stresses  at  t=200  sec.  ranging  from  30  MPa  (4.4  ksi)  to  80  MPa  (1 1 .6  ksi)  could 
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well  exceed  a typical  tensile  strength  of  high-strength  concrete.  Failure  to  include  the 
effects  of  moisture  in  the  process  of  determining  temperature  distribution  results  in 
numerical  predictions  of  thermal  gradient  stresses  that  are  nearly  fifty  per  cent  greater 
at  early  stages  of  thermal  loading  (refer  to  Chapter  5).  This  over-estimation  of  thermal 
gradient  stress  values  becomes  more  exaggerated  as  thermal  loading  continues. 
Therefore,  the  influence  of  moisture  migration  and  evaporation  on  heat  transfer  must 
be  considered  in  order  to  obtain  valid  numerical  predictions  of  thermal-dilatation 
stresses. 

The  effects  of  temperature  distribution  combined  with  boundary  restraints  on 
hydro-thermo-elastic  stresses  have  been  identified.  Thermal  gradient  stresses  increase 
with  temperature  gradient.  However,  the  temperature  prediction  by  a conventional 
transient  thermal  analysis  is  higher  than  by  the  coupled  moisture  and  heat  transport 
analysis.  This  leads  to  even  steeper  development  of  thermal  gradients  and  thus  causes 
overestimation  of  thermal  gradient  stresses.  Numerically  predicted  and 
macroscopically  volume-averaged  pore  pressure  has  nearly  negligible  effects  on 
development  of  the  maximum  principal  stress  during  the  first  five  minutes  of  fire 
exposure.  Although  magnitudes  of  pore  pressure  buildup  are  much  lower  that  those  of 
thermal  gradient  stresses,  one  should  not  simply  neglect  the  effects  of  moisture  in 
thermal  analysis.  It  has  been  found  that  moisture  movement  affects  temperature  rise 
and  gradient  development.  Therefore,  coupled  heat  and  mass  transfer  processes  are 
involved  in  the  development  of  pore  pressure  establish  the  thermodynamic  states  of 
high-strength  concrete  when  exposed  to  elevated  temperatures  (Chung  and 
Consolazio  2004). 
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It  is  noted  that  the  magnitudes  of  thermo-elastic  gradient  stresses,  predicted  by 
the  present  finite  element  stress  analysis  model,  in  the  near  surface  of  the  columns 
may  be  overestimated.  Young’s  modulus  of  elasticity  of  cover  concrete  decreases  as 
temperature  increases  (Phan  et  al.  2001,  Janotka  and  Bagel  2002).  The  combined 
effect  of  thermal  strains  that  increase  as  temperature  differentials  increase  and 
Young’s  modulus  of  elasticity  of  concrete  that  decreases  at  an  elevated  temperature 
needs  to  be  investigated  for  more  accurate  assessment  of  thermal  dilatation  induced 
stresses.  Reduction  of  concrete  strength  in  both  compression  and  tension  should  be 
considered  in  order  to  determine  the  locations  and  time  of  spalling.  The  assessment  of 
compatibility  needs  to  be  investigated  further  in  regard  to  1 ) constituents  of 
reinforced  concrete  and  2)  the  geometry  of  the  structural  members.  For  example,  the 
differences  in  thermal  expansion  between  cement  paste  and  aggregate,  and  between 
mortar  and  steel  reinforcement  could  cause  interfacial  bond  failure.  For  a rectangular 
column,  thermal  expansion  will  be  smaller  along  a shorter  dimension  of  the  cross- 
section  than  along  the  longer  dimension  of  the  cross-section.  Thus,  an  increase  in 
compression  in  the  longer  dimension  will  result  whereas  a decrease  in  compression  or 
even  tensile  forces  can  result  in  the  shorter  direction  (Usmani  et  al.  2001). 

6.2  Contribution  of  Pore  Pressure  at  a Later  Stage  of  Heating 

Another  scenario  of  spalling  can  be  drawn  when  pore  pressures  become  quite 
substantial  under  a sustained  thermal  loading  condition.  Potential  influence  of  pore 
pressures  on  spalling  has  been  investigated  under  an  assumption  that  the  concrete 
cover  remains  in  place  until  the  pore  pressures  build  up  significantly. 
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Figures  6-1  and  6-2  illustrate  the  thermodynamic  field  variables  developed  in 
two  different  locations  after  time  t=600sec.  (10  minutes)  when  all  the  four  sides  of  the 
square  column  discussed  in  Chapter  5 (refer  to  Figure  5-6)  are  initially  exposed  to  a 
hydrocarbon  pool  fire  modeled  using  ASTM  El  529.  In  contrast  to  negligible  values 
of  pore  pressure  at  an  early  stage  of  heating  (i.e.,  t=200  sec.),  numerically  computed 
volume-averaged  pore-pressure  buildup  in  a plane  where  shear  reinforcement  is 
located  (Figure  6-la)  exceeds  4 MPa  (0.58  ksi)  at  time  t=400  sec.  and  reaches  5.5 
MPa  (0.8  ksi)  at  time  t=600  sec.  (Figure  6-lb).  Interestingly,  the  development  of  peak 
pore  pressure  is  predicted  along  the  interfaces  between  concrete  and  shear 
reinforcement  parallel  to  the  surfaces  exposed  to  fire. 

Similar  trends  of  pore  pressure  are  found  in  a plane  located  at  20  mm  vertically 
above  the  shear  reinforcement  (Figure  6-2a).  Occurring  in  the  comer  region  of  the 
square  column,  the  maximum  value  of  volume-averaged  pore  pressure  in  Plane  B 
exceeds  6.7  MPa  (0.97  ksi)  at  time  t=600sec.  (10  minutes). 

Unexpectedly,  as  heating  continues,  higher  pore  pressure  is  developed  in  the 
concrete  region  between  adjacent  shear  ties  (Plane  B)  in  the  column  than  in  the 
vicinity  of  the  shear  reinforcement  (Plane  A).  The  shear  reinforcement  was  indeed 
blocking  horizontal  moisture  flow  in  the  xz-plane  and  forced  the  moisture  to  permeate 
vertically  in  the  directions  parallel  to  the  column  surfaces.  According  to  experimental 
works  by  Elvery  et  al.  (1973)  and  Butler  (1981),  the  pore  pressure,  averaged  using  the 
effective  porosity  at  tensile  failure  of  high-permeability  concrete,  was  measured  in  a 
range  of  1 .5-2.5  MPa.  The  pore  pressure  alone  predicted  by  the  present  study  could 


cause  tensile  failure  of  concrete. 
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(a)  Finite  difference  model 
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Figure  6-1  Thermodynamic  states  of  a square  column  segment  at  t=600  sec.  (data 

plotted  for  Plane- A of  Figure  5-6) 


In  addition,  the  pore  pressure  distributions  with  respect  to  formation  of 
saturation  fronts  are  identified  in  Figures  6- Id  and  6-2d.  It  is  noted  that  locations  of 
maximum  pore  pressure  buildup  are  found  to  be  within  a two-phase  zone  of  finite 
length  away  from  the  column  surface  adjacent  to  locations  where  the  degree  of 
saturation  maximizes  (Pruess  et  al.  1987).  Very  similar  temperature  distributions  are 
observed  in  both  Plane  A and  Plane  B (Figures  6-lb  and  6-2b). 
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(a)  Plane-B  mesh  layout 


(b)  Pore  pressure 
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Figure  6-2  Thermodynamic  states  of  a square  column  segment  at  t=600  sec.  (data 

plotted  for  Plane-B  of  Figure  5-6) 


Comparing  Figure  6- Id  to  Figure  6-2d,  we  observe  more  prominent 
development  of  moisture  clogging  near  the  shear  reinforcement  in  Plane  A.  Note  that 
this  liquid  water  saturation  is  just  barely  below  the  maximum  theoretical  value  of 
unity.  A wall  of  moisture  constituting  a form  of  “heat  sink”  is  formed  as  shown  in 
Figure  6-3.  Therefore,  a large  amount  of  thermal  energy  is  consumed  in  the 
evaporation  process  of  liquid  water  that  is  retained  by  moisture  clogging  and 
subsequently  temperature  rise  in  the  moisture-clogged  region  is  retarded.  This  results 
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in  a significant  reduction  of  thermal  gradients  and  thus  magnitudes  of  thermal 
gradient  stresses  whereas  pore  pressure  continues  to  increase.  It  is  concluded  that  at  a 
later  stage  of  heating  when  moisture  clogging  becomes  prominent,  the  contribution  of 
pore  pressure  to  spalling  is  significant.  Violent  fracture  failure  (i.e.,  explosive 
spalling)  of  high  strength  concrete  has  been  reported  during  the  first  10-20  minutes 
of  building  fire  exposure  in  field  and  laboratory  tests  (Kodur  1 999,  Chana  and  Price 
2003).  In  their  tests,  large  areas  of  steel  reinforcement  were  exposed  and  concrete 
continued  to  spall.  With  decreasing  severity,  this  explosive  spalling  was  observed 
until  25  minutes  after  ignition.  The  combined  effects  of  thermal  dilatation  induced 


stresses  and  pore  pressure  result  in  a continuation  of  explosive  spalling. 


0 o 


Figure  6-3  Moisture  clog  formed  in  a square  column 


One  factor  influencing  pore  pressure,  which  has  not  been  investigated  here,  is 


the  significant  increase  of  material  flow  parameters  such  as  the  permeability  of 
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concrete  at  high  temperatures  (Phan  and  Carino  1998,  Janotka  and  Bagel  2002). 
Expansion  of  micropore  volume  and  cracking  induced  by  shrinkage  (occurring  on 
loss  of  moisture  due  to  dehydration  of  cement  gels  at  high  temperatures)  are 
phenomena  that  cause  a sudden  increase  of  the  permeability  of  concrete.  Thus,  the 
values  of  the  pore  pressure  predicted  by  present  models  may  be  higher  than  those 
actually  occurring  in  high-strength  concrete  (Consolazio  et  al.  1998,  Kalifa  et  al. 
2000).  Nonetheless,  the  material  model  developed  in  this  study  can  be  further 
improved  when  additional  experimental  data,  particularly  permeability  of  high- 
strength  concrete  at  elevated  temperatures,  are  available.  Completely  defining  the  role 
of  pore  pressure  in  initiating  spalling  will  require  additional  investigation  and 
experimental  testing. 


CHAPTER  7 
CONCLUSION 


An  improved  quantitative  understanding  of  multi-dimensional  states  of  stress  inside 
reinforced  concrete  structural  elements  under  severe  fire  condition  has  been  obtained  in 
this  study.  The  combined  finite  difference  and  finite  element  model  developed  in  this 
study  provided  valuable  insights  into  hydrothermal  behavior  of  reinforced  high-strength 
concrete  structural  elements.  By  developing  thermo-elastic  stress  analysis  techniques  for 
concrete  structural  elements  exposed  to  severe  fire,  this  research  has  focused  on  the 
development  of  a practical  computer  methodology  that  can  be  used  to  study  concrete 
thermal  spalling.  Thermal  mechanisms  such  as  thermal  expansion  and  thermal  bowing 
have  been  found  to  produce  severe  stress  states  (large  compressive  or  tensile  stresses) 
depending  on  the  temperature  distribution  in  the  concrete  member  and  boundary  restraint 
conditions. 

A comprehensive  theoretical  understanding  of  how  elevated  temperatures  generated 
during  a fire  affect  the  thermodynamic  state  of  concrete  has  been  established  during  this 
research.  By  establishing  practical  models  that  take  into  account  the  fluid  flow  and 
thermal  properties  of  concrete,  the  effects  of  coupled  moisture  and  heat  transport 
phenomena  occurring  in  heated  high-strength  concrete  have  been  studied. 

Components  of  the  research  involved  studies  of  material  characteristics  of  concrete 
as  a porous  medium  and  effects  of  multiphase  flow  on  magnitude  of  pore-pressure  build 
up,  development  of  thermal  gradients,  and  evolution  of  degree  of  saturation.  Also,  the 
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effects  of  thermal  exposure  conditions  such  as  rate  of  heating  and  thermodynamic 
conditions  of  concrete  in  ambient  conditions  (e.g.,  initial  degrees  of  saturation  at  various 
levels  of  relative  humidity)  were  examined. 

Using  thermodynamic  state  data  obtained  from  the  computer  finite  difference 
models,  the  influence  of  thermally  induced  stresses  on  spalling  has  been  investigated  in 
consideration  of  the  location  of  steel  reinforcement,  geometrical  structural  shapes,  and 
constraints  on  the  structural  element.  In  the  following,  concluding  remarks  of  this  study 
are  stated  in  regard  to  the  influence  of  thermally  induced  stress  on  spalling.  Discussion  of 
future  research  areas  in  relation  to  the  numerical  methodology  proposed  in  this  study  is 
then  provided. 


7.1  Thermal  Spalling  of  Concrete 

The  effects  of  thermal  exposure  conditions  such  as  rate  of  heating  and  initial 
thermodynamic  conditions  of  concrete  (e.g.,  saturation)  have  been  found  to  be  significant 
factors  that  influence  moisture  flow  inside  heated  concrete.  The  important  role  of 
moisture  has  been  identified  by  simulation  results  of  thermal-dilatation  stresses  in 
reinforced  concrete  columns.  Due  to  low  conductivity  and  high  heat  capacity  of  partially 
saturated  concrete,  the  rate  of  heat  transfer  is  very  low.  Accumulation  of  thermal  energy 
results  in  the  development  of  a steep  thermal  gradient  in  the  concrete  cover,  which  can  be 
even  more  intensified  by  rapid  rate  heating.  A large  amount  of  thermal  energy  is 
consumed  by  the  evaporation  process  of  moisture  in  concrete  pores.  Maximum 
temperature  rise  in  the  near-surface  region  is  affected  as  is  the  development  of  the 
associated  temperature  gradient. 
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As  detailed  in  earlier  chapters,  after  a few  minutes  of  fire  exposure,  moisture 
clogging  occurs.  Severe  and  early  formation  of  moisture  clogging  is  found  in  the  case  of 
high-strength  concrete  more  than  in  the  case  of  normal  strength  concrete  exposed  to  fire. 
Several  factors  including  lower  permeability  and  higher  degree  of  initial  saturation  at  an 
ambient  condition  create  a pseudo  heat  sink  and  impermeable  boundary  inside  high- 
strength  concrete.  Thus,  the  formation  of  a moisture  clog  is  closely  related  with,  and 
directly  affects  the  rate  of  temperature  rise  and  pore-pressure  development. 

The  influence  of  pore-pressure  on  thermal  spalling  of  high-strength  concrete  has 
been  examined.  By  means  of  linear  thermoelastic  analysis,  the  contribution  of  the  pore 
pressure  buildup  to  the  effective  stresses  during  the  first  five  minutes  of  fire  exposure  has 
been  found  to  be  negligible  compared  to  that  of  thermal  dilatation  stresses.  Based  on 
stress  analysis  results  presented  herein,  thermally  induced  spalling  appears  to  be 
primarily  driven  by  thermal  dilatation  stresses,  not  by  pore  pressure.  Development  of 
thermal  gradient  stresses  has  been  found  to  be  significantly  influenced  by  boundary 
restraint  conditions.  Another  important  factor  that  could  affect  thermal  spalling  of 
structural  concrete  members  has  also  been  investigated.  Restraint  conditions  have  been 
found  to  have  a significant  effect  on  the  development  of  internal  forces  and  the 
displacements.  For  instance,  boundary  constraints  such  as  vertical  restraints  imposed  on 
the  3-D  column  segment  model  and  the  2-D  axisymmetric  column  model  produce  large 
tension  forces  when  steep  thermal  gradient  is  developed  in  the  outer  hotter  region  of  the 
columns.  For  the  source  of  structural  restraints,  it  is  not  so  clear  whether  sufficient 
restraints  would  remain  available  as  temperature  rises.  This  is  an  area  that  needs  to  be 
investigated  further. 


174 


Effects  of  moisture  on  development  of  thermal  stresses  have  been  confirmed  and 
considered  in  multi-dimensional  finite  element  stress  analyses.  Thus,  the  effects  of 
moisture  movement  on  thermal  stresses  should  be  considered  in  determination  of  the 
effective  stresses  due  to  1)  influence  in  pore  pressure,  2)  potential  association  with  crack 
growth  in  the  fracture  processes,  and  most  importantly  3)  determining  the  distribution  of 
temperature. 

It  must  be  noted  that  within  this  study,  any  prediction  of  material  failure  involved  in 
thermal  spalling  would  be  limited  by  the  present  linear  elastic  stress  analysis. 
Temperature  dependency  of  the  key  material  thermo-mechanical  parameters  and  failure 
mechanisms  associated  with  the  material  strength  parameters  should  be  considered  in 
future  research  to  quantitatively  determine  more  accurate  stress  fields.  In  addition,  in 
order  to  determine  whether  a structural  element  will  be  safe  under  a severe  thermal 
loading,  the  state  of  stress  should  be  computed  at  all  critical  points  of  the  structural 
element,  i.e.,  at  all  points  where  stress  concentrations  are  likely  to  occur.  This  may  be 
done  in  a number  of  ways  by  using  stress  concentration  factors  when  the  theory  of 
elasticity  is  used  to  determine  the  state  of  stress  at  each  critical  point. 

When  material  failure  induced  by  macroscopic  cracks  is  considered  in  a structural 
element,  special  care  needs  to  be  taken.  When  a crack  is  initiated  due  to  an  initial  fracture 
in  a thermally  loaded  structural  element,  it  is  necessary  to  determine  whether  the 
macroscopic  crack  will  tend  to  propagate  under  the  expected  loading  condition  and  cause 
structural  failure,  or  whether  it  will  remain  stable.  This  requires  a computational  analysis 
involving  the  fracture  energy  associated  with  the  growth  of  the  crack  and  inclusion  of 
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both  stress  and  internal  pore  pressure  effects.  Such  an  analysis  is  beyond  the  scope  of  this 
text  and  should  be  carried  out  by  the  methods  of  fracture  mechanics. 

7.2  Future  Research  Work 

This  research  has  considered  key  issues  related  to  fire  induced  spalling  of 
reinforced  concrete  structural  elements  in  extreme  fire  events.  A new  numerical 
methodology  is  presented  to  analyze  heat  and  mass  transport  phenomena  in  concrete  at 
high  temperatures  that  results  in  predictions  of  mechanical  behavior.  The  computational 
steps  associated  with  the  development  of  the  numerical  models  can  be  further  improved 
by  considering 

1 . Computer  modeling  of  changes  in  mechanical  properties  of  high-strength 
concrete  such  as  modulus  of  elasticity,  tensile  and  compressive  strength, 
Poisson’s  ratio,  linear  thermal  expansion  coefficient,  and  heat  capacity 
coefficient  at  elevated  temperatures,  should  be  included  in  future  models. 

2.  Computer  modeling  of  changes  in  thermophysical  properties  such  as 
permeability,  thermal  conductivity,  specific  heat,  and  density  at  high 
temperatures  should  be  considered  in  the  transient  thermal  analysis. 

3.  The  model  can  be  extended  to  predict  material  failure  at  high  temperatures 
so  as  to  determine  the  fire  resistance  of  high-strength  concrete.  The  effects 
of  expansion  of  aggregates  and  shrinkage  of  cement  paste  can  be 
incorporated  in  the  thermo-mechanical  stress  analysis. 

4.  The  effects  of  crack  opening  by  thermal  gradient  stresses  should  be 
considered  with  the  theory  of  thermofracture  mechanics.  Genuine  efforts 
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must  be  made  to  include  the  influence  of  the  pore-pressure  buildup  on 
tension  cracking  in  the  analytical  solution  for  the  governing  equations. 

The  integrated  computer  analysis  tool  developed  in  this  study — integration  of 
TOUGH2  and  PLASFEM — will  serve  as  a basis  for  building  an  analytical  model  capable 
of  modeling  progressive  material  failure  that  has  been  found  to  be  a significant 
component  of  spalling  phenomenon  in  concrete  structures  exposed  to  fire. 
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